NASA TECHNICAL NOTE 


C s 4 



ar;;„ 

KIRTLAND Ai 


un 


NUMERICAL CALCULATIONS FOR THE 
CHARACTERISTICS OF A GAS FLOWING 
AXIALLY THROUGH A CONSTRICTED ARC 


by Velvin R. Watson and Eva B. Pegot 

Ames Research Center 
Moffett Field, Calif. 



\?2] 


s?y 






<-c 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHINGTON, D. C. • 


I 


1 1 


D-40 42 






JUNE 1967 

t 


TECH LIBRARY KAFB, NM 



TECH LIBRARY KAFB, NM 


H 


□ 13D77 C 1 

NASA TN D-4042 


NUMERICAL CALCULATIONS FOR THE CHARACTERISTICS OF A GAS 

FLOWING AXIALLY THROUGH A CONSTRICTED ARC 

By Velvin R. Watson and Eva B. Pegot 

Ames Research Center 
Moffett Field, Calif. 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


For sale by the Clearinghouse for Federal Scientific and Technical Information 
Springfield, Virginia 22151 — CFSTI price $3.00 



NUMERICAL CALCULATIONS FOR THE CHARACTERISTICS OF A GAS 


FLOWING AXIALLY THROUGH A CONSTRICTED ARC 

By Velvin R. Watson and Eva B. Pegot 
Ames Research Center 


SUMMARY 


Numerical programs to obtain solutions for the characteristics of a gas 
flowing axially through a constricted axe are presented. The numerical pro- 
grams use real equilibrium gas properties and solve simultaneously the energy, 
momentum, and continuity equations. Axial conduction, radial pressure gradi- 
ents, and radial voltage gradients are neglected. The solutions give the arc 
characteristics in sufficient detail to evaluate the many approximate solu- 
tions, and the computing time (approximately 2 min) is sufficiently short 
that the programs may be used directly to obtain design criteria for plasma 
generators. 

Numerical solutions for arcs within O.635 and 1 . 27 -cm- diameter 
constrictors are presented and the design optimization of constricted-arc 
plasma generators is discussed. The numerical solutions indicate that with 
nitrogen, total enthalpies in excess of 5x10 s j/kg and velocities in excess 
of 18,000 m/s may be obtainable at the exit of a constricted-arc plasma 
generator. 


INTRODUCTION 


The constricted arc has recently been employed (as in the device shown 
in fig. 1) to generate hot, dense plasma flows for production of the very 
high heat fluxes required in materials testing, and to produce thrust. 
Approximate solutions for the characteristics of the flow through this con- 
stricted arc (refs. 1-6) predict the arc -column characteristics with varying 
degrees of accuracy. The discrepancies between actual and predicted charac- 
teristics are a result of the simplifying approximations incorporated into 
each model to permit an analytical or semianalytical solution; the extent to 
which these various approximations influence the behavior of the solutions, 
however, is not readily apparent. 

The purpose of this paper is to present numerical methods for solving a 
more complete model with fewer and, presumably, more realistic assumptions. 
These solutions are sufficiently detailed to evaluate the simplified models 
and to gain further insight into the behavior of the constricted arc. 

An evaluation of the approximations for the theoretical model presented 
in reference 1 has already been made and was presented in reference 7- The 
numerical solutions used for this evaluation are presented in appendix A. 



Several arc characteristics were studied with these numerical methods 
and the results of these studies are presented. The results indicate methods 
for optimizing the design of a plasma generator for particular arc character- 
istics, and these methods are discussed. 

In a previous work by Masser (ref. 8) a similarly complete method was 
presented for solving the constricted -arc problem. However, it will be shown 
that to obtain the same accuracy Masser f s method may require a longer comput- 
ing time than the method presented herein. 


NOMSNCLA.TUHE 

A cross-sectional area of the constrictor, m 2 

A- parameter of the linear approximation a = A <p, mho/W 

O G> 

a zero frequency speed of sound, m/s 

c specific heat at constant pressure, J/kg °K 

Jr 

d diameter 

E voltage gradient , V/m 

E$o 2.4/RA g 1/2 

mass -average enthalpy, J/kg (H m = 0 at 0° K) 

H g enthalpy averaged over space, J/kg (H s = 0 at 0° K) 

theoretical mass average enthalpy in the asymptotic region of the arc 
(ref. 1), 0.133 c p l/kRA g 1/2 , J/kg 

h enthalpy, J/kg (h = 0 at 0° K) 

hp center -line enthalpy, j/kg 

iLj. total enthalpy, J/kg 

hoo theoretical center -line enthalpy in the asymptotic region of the arc 
(ref. 1), 0.307 epl/kRA g 1/2 , j/kg 

I current, A 

k thermal conductivity, W/m °K 

tii mass -flow rate, kg/s 

p pressure, N/m 2 


2 


p Q constrictor inlet pressure, N/m 2 

q local heat transfer rate from the surface of the arc column to the 
constrictor wall, W/m 2 

q r radiative component of the heat transfer rate to the constrictor wall. 
W/m 2 

q-fc total heat transfer rate to the constrictor wall, w/m 2 

qco theoretical heat transfer rate from the surface of the arc column in the 
asymptotic region of the arc (ref. 1 ), O.383 l/R 2 A g l/2 , w/m 2 

R constrictor radius, m; also gas constant 

r radial distance from the axis of the arc column, m 
o 

T temperature, K 

u axial velocity, m/s 

v radial velocity, m/s 

Z compressibility 

z axial distance along the column, m 

z 0 characteristic length for the arc (ref. 1 ), mCp/:tk, m 
0 azimuthal coordinate of the arc column, radians 
\i viscosity, Ns/m 2 

P 0 magnetic permeability, n/A 2 

p density, kg/m 3 

a electric conductivity, mho/m 

cp thermal conductivity function, J k dT, W/m (cp = 0 at 0° K) 

DEVELOPMENT OF NUMERICAL ANALYSIS 


Theoretical Models for the Numerical Calculations 

The theoretical models for this work are assumed to be governed by one 
of the following two sets of equations. Equation (l) is applicable to the 
case of axial symmetry; equation (2) allows for asymmetric flow without 
swirl. 
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In presenting the equations in the above form it is assumed that: 

1. The gas flow is steady and laminar. 

2. The electric discharge is stationary and the electric potential is 
constant on planes perpendicular to the axis. 

3- Axial heat conduction is negligible compared to radial heat 
conduction. 

4. Lorentz forces are negligible compared to dynamic and static 
pressure forces. 

5- The radial pressure gradient is negligible compared to the static 
pressure. 

A detailed discussion of these assumptions is presented in reference 1. 

In comparison, the model of reference 1 contains all of the above 
assumptions and approximations in addition to the following simplifying 
approximations : 
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1- The mass flux is assumed constant throughout the constrictor. 

2. The enthalpy profile at the constrictor inlet is assumed to be a 
Bessel function. 

3- The more important gas properties - enthalpy, thermal conductivity 
potential, and electrical conductivity - are assumed to be linearly related 
and the radiance and viscosity are neglected. 

The gas properties used in the numerical solutions for this report are 
theoretical estimates of real, equilibrium properties for nitrogen and hydro- 
gen and are shown in figures 2 and 3- (in some of the numerical calculations 
the radiance of air was used in place of the radiance of nitrogen; so the 
radiance of air is also shown in figs. 2 and 3*) Table I lists the source 
for each of these gas properties (refs. 9-H) • 


Numerical Procedures 

The numerical procedure for solving equations (l) and (2) is to satisfy 
finite-difference representations of the equations by forward marching from 
assumed upstream and constrictor wall boundary conditions. The step-by-step 
procedure for obtaining solutions to equation (2) is presented below. Fig- 
ure k illustrates the arrangement of the finite-difference network. 

1. Establish the following initial and boundary conditions: 

(a) Enthalpy distribution at station 1 

(b) Velocity distribution at station 1 

(c) Enthalpy of the gas adjacent to the constrictor wall, that is, 

the enthalpy of the gas at the wall temperature 

(d) The pressure at station 1 

(e) The total current through the constrictor 

2. Evaluate the distribution of the following gas properties at sta- 
tion 1 from the known values of enthalpy and pressure: 

(a) Thermal conductivity potential 

(b) Electrical conductivity 

(c) Radiance 

(d) Viscosity 

(e) Density 
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3 • Evaluate the electrical conductance at station 1 by integrating the 
electrical conductivity over the area at station 1 . Evaluate the electrical 
voltage gradient at station 1 from the total current and electrical 
conductance. 

4. Evaluate the total flow rate by integrating the local mass flow, pu, 
over the area at station 1 . 

5* Compute the enthalpy at all mesh points at station 2 from the energy 
equation. 

6 . Estimate a trial pressure drop between stations 1 and 2 ; calculate 
the velocity at station 2 from the momentum equation* 

7- Evaluate the density for all mesh points at station 2 from the known 
values of enthalpy and pressure. 

8 . Evaluate the total flow rate at station 2 by the same method as 
described in step 4. 

9. If the total flow rate at station 2 agrees with the specified flow 
rate to within the specified accuracy, continue by repeating steps 2 through 
9 for each new station; otherwise, change the trial pressure drop between 
stations 1 and 2 and repeat steps 6 , 7 , and 8 . (To increase the flow rate at 
station 2 , the trial pressure drop must be increased for subsonic flow values 
and decreased for supersonic flow values.) 

If a repeated change of the trial pressure drop does not produce a 
calculated flow rate sufficiently near the specified flow rate, then either 
the calculations have become unstable or the flow has been aer ©dynamically 
choked in the tube. If the flow has become choked, the calculations can be 
continued into the supersonic region by causing the diameter to increase and 
by reversing the conditions on the trial pressure at the station of choking; 
that is, if the flow rate is too large at the station, then the trial pressure 
drop should be increased. 

The procedure for solving equation (l) is similar except that the radial 
convection term is present in the momentum and energy equations. The radial 
mass flux required for the radial convection terms is determined by means of 
a local mass flux balance in each volume element, starting with the center 
element; that is, for the center volume element between stations n and n + 1 , 
the radial mass flux into this volume is equal to the axial mass flux out at 
station n + 1 minus the axial mass flux in at station n. Eor computing 
convenience the radial momentum and energy fluxes between stations n and 
n + 1 are added downstream between stations n + 1 and n + 2 . 


Fortran Programs 

The programs that solve equations (l) and (2) are written in Fortran II 
and are presented in appendixes B and C. Included in the appendixes are 
descriptions of the required subroutines and definitions of the variables used 
within the programs . 
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The programs that calculate the gas properties from known values of 
enthalpy and pressure are given in appendix D. 

The program that prepares the magnetic tape required for the program 
given in appendix D from a table of gas properties is given in appendix E. 


Display of Solutions 

The above Fortran programs solve for the local state properties and 
velocity of the gas throughout the constricted arc. For the symmetric arc, 
the distribution of these local properties can easily be visualized with an 
oblique projection as shown in figure 5- Figure 5(a) represents a constant 
area constrictor whereas figure 5(b) represents a constrictor with a contoured 
inlet section. The horizontal scale represents the radial distance from the 
axis of symmetry, the oblique scale represents the axial distance within the 
constrictor, and the vertical scale represents the magnitude of the local 
property - in this illustration, enthalpy. The local values of enthalpy, 
mass flow, energy flow, velocity, and momentum flow are illustrated with these 
oblique projections in parts (a) through (e) of figures 6 through 3^** 

All arc -column characteristics of interest can be calculated from the 
local properties obtained from the numerical solutions* Local properties 
that have been calculated and plotted as functions of axial distance are 
center -line enthalpy, average enthalpy (averaged over both space coordinates 
or mass flux coordinates) voltage, constrictor wall heat transfer rate, and 
static pressure. These are shown in parts (f) through (k) , respectively, of 
figures 6 to 3^* The predictions from the analytical model of reference 1 
for the constricted arc with negligible radiation losses are also shown in 
parts (f) through (j) for comparison. The program that plots these graphs 
from the results of the numerical program given in appendix B, is given in 
appendix F. 

The results for the asymmetric arc cannot be shown completely by the 
oblique projections used to display the results for the symmetric arc. 
Nevertheless, the tabular output of the above programs can be visualized if 
the results are arranged as shown in figure 35> wherein the magnitudes of the 
enthalpy and the mass flux are shown as functions of radius and azimuthal 
position at each axial station. 

Although the arc is not symmetric, the hot section of the arc does not 
rotate with axial distance, and the values of the gas properties on a plane 
of constant azimuthal position that passes through the hot spot is of inter- 
est. The values of the gas properties on this plane can be shown by the 
oblique projection as illustrated in figure 3 6 . 


Computational Accuracy 

Mesh size . - The calculations for the symmetric arc were made with 13 
radial mesh points between the constrictor axis and the constrictor wall and 
with sufficient axial stations to yield a stable solution (usually between 
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200 and 1000 stations) . A comparison was made for a case in which only the 
mesh size was changed. The number of radial mesh points was changed from 13 
to 26 and the number of axial stations was changed from 263 to 8^7 • The heat 
flux to the constrictor wall, which is obtained from the derivatives of the 
local enthalpy profile, changed the most (nearly 10 percent) whereas the mass- 
averaged enthalpy which is obtained from the integral of the local enthalpy 
profile, changed by less than 3 percent. The heat transfer rate is the only 
property that is obtained from a derivative, and most of the arc properties 
changed less than 5 percent as a result of the change in mesh size. 

Stability . - A rigorous stability criterion for the finite difference 
solutions of these nonlinear equations is not known. However, a stability 
criterion for the finite difference solution of linear parabolic equations is 
that 


Az ua 
£r 2 2 


where 

u gas velocity 
a thermal diffusivity 


Furthermore, whenever the local values of Az/Ar 2 exceeded this limit in the 
numerical calculations of the nonlinear equations, instabilities were encoun- 
tered. Therefore, the above criterion for stability appears to apply locally 
for the nonlinear equations, and a step was added to the above numerical pro- 
cedure to keep the local value of Az/Ar 2 less than this limit. 

Instabilities were encountered in the program for the numerical calcula- 
tions for equation (l) when the pressures were high. These instabilities were 
eliminated by changing the method of handling the radial convection. The 
instabilities were encountered when the momentum and energy that were con- 
verted radially between stations n and n + 1 were added at station n + 2 as 
described under numerical procedures. Adding one -fourth of this radially con- 
vected momentum and energy to each of the stations n+2,n+3j n+ ^j and 
n + 5 eliminated the instability. 

Computation time .- The programs presented herein are written in Fortran 
II and were executed by an IBM 7094 . The computation time required for the 
calculation of the solutions is approximately 2 minutes. The time for each 
case shown in figures is given in table II. 

The computing time is proportional to the dimensionless length of the 
constrictor and is approximately proportional to the cube of the number of 
radial mesh points. 
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Discussion 


It will be shown in the sections that follow that the numerical solutions 
display the qualitative features of the constricted arc in sufficient detail 
to gain an understanding of the arc column and to evaluate the various 
approximate analyses of the constricted arc. 

The computation times for the programs for the first set of equations 
(approximately 2 min) are sufficiently short that it appears feasible to use 
these programs directly to obtain design criteria for plasma generators. 

Several physical phenomena that were neglected in these calculations can 
be included without greatly increasing the complexity or calculation times 
of the program. For example, if the operating conditions are changed so that 
the magnetic pressure term will not be negligible, this term can easily be 
included in the calculations . (The current distribution is known and because 
of symmetry, the magnetic field is easily calculated.) Also, the effects of 
radiation absorption could be approximated if the frequency of the radiation 
were divided into intervals and the radiation and absorption path length were 
specified for each frequency interval. The radiation at all frequencies for 
which the path length is much greater than the constrictor radius could be 
specified with one term and the gas considered to be transparent to this 
radiation, and the radiation at all frequencies for which the path length is 
much less than the radial mesh increment could be included as a single term 
similar to the thermal conduction term. Since the variation of radiation 
with axial distance is small compared to the variation of radiation with 
radial distance, the radiation absorbed at each axial station can be calcu- 
lated approximately from the temperature and density profiles of the previous 
axial station. For the higher flow rates, turbulence may become significant 
and the convective transfer of energy and momentum could be approximated with 
the Prandtl mixing length theory. 

The starting conditions (i.e., the enthalpy and velocity at the first 
axial station) are usually unknown and must be assumed. Figures 6 through 9 
compare solutions in which only these assumed starting conditions were 
changed. Here one can see that the effect of the starting conditions on the 
solutions at very short characteristic lengths (i.e., very small z/z G ) is 
severe. Nevertheless, this effect diminishes rapidly downstream, and even at 
moderately small characteristic lengths the effect of the starting conditions 
is small. 

A procedure for the numerical solution of < equation (l) is given by 
Masser (ref. 8). Comparing Masser's work with the present procedure, the two 
differ in their respective location of the radial mesh points; Masser 1 s are 
fixed on streamlines whereas in the present program they are fixed in space. 
Since the dominant heat -transfer mechanisms (i.e., thermal conduction, radia- 
tion, and ohmic heating) are dependent upon space coordinates rather than 
mass flux coordinates, the accuracy with which these terms can be calculated 
depends upon the mesh distances in space coordinates. From part (b) of the 
figures illustrating the solutions, it can be seen that the mesh distances in 
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space would vary radically if the mesh points were fixed on streamlines. 
Masser T s procedure, therefore, would probably require more mesh points for 
the same degree of accuracy. 


EVALUATION OF THE NUMERICAL SOLUTIONS 


The numerical solutions were compared with an exact analytical solution 
to determine whether the numerical procedures yield the proper solutions to 
the equations for the theoretical models. Then the numerical solutions were 
compared with experimental measurements to determine whether the theoretical 
model yields the correct arc column features. 


Comparison of the Numerical Solution With an Exact Solution 

An exact analytical solution for a theoretical model of the arc column 
with one set of boundary conditions is given in reference 1. The theoretical 
model of reference 1 is the same as the theoretical model used for the numer- 
ical calculations for this paper except that the model of reference 1 contains 
additional simplifying approximations. (The list of additional simplifying 
approximations Is given in a previous section "Theoretical Models for the 
Numerical Calculations.") In order to compare the numerical solutions with 
the exact solution, the additional approximations of reference 1 were incor- 
porated into the model used for the numerical calculations (e.g., linearized 
gas properties were used for the numerical calculations) . The number of 
radial mesh points for the numerical calculation was 13 and the number of 
axial mesh points was 195- The difference between the numerical solution and 
the exact solution was less than 1 percent (nearly Indistinguishable on a 
graphical illustration) indicating that the numerical procedures give the 
proper solution to the equations for the theoretical model. 


Comparison of the Numerical Solutions With Experimental Results 

The numerical calculations are compared with the experimental measure- 
ments of the mass average enthalpy, voltage gradient, and wall heat transfer 
rate in a 1.27 -cm-diameter constricted -arc plasma generator (figs. 10-14), 
and in a 0. 63 5 -cm-diameter constricted-arc plasma generator (figs. 15-20). 
(Figs. 10-20 display 10 arc column features even though there are experimental 
measurements for only three of these features. The three features for which 
the numerical solutions are compared with experimental measurements are shown 
in parts (h) , (i) , and (j) of these figures.) 

The mass average enthalpy was measured experimentally by subtracting the 
heat losses from the electrical power Input. The 1.27 -cm-diameter constrictor 
was made in two modules and only the average heat transfer rate to each module 
was measured. Therefore, the mass average enthalpy could be determined only 
at the end of each module. The electrical voltage gradient was not measured 
for all runs. 

The starting conditions for the numerical calculations are not known 
a priori; the distribution of enthalpy and velocity at the inlet must be 
assumed. For the 1.27 -cm-diameter constrictor, an arbitrary and fairly low 
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value of enthalpy was selected at the start. When the average enthalpy in 
the numerical calculations became equal to the average enthalpy of the gas 
entering the constant area constrictor of the plasma generator, the axial 
distance of the numerical calculations was made to coincide with the start of 
the constrictor in the plasma generator. This is shown graphically in part 
(h) of the figures wherein the curve for the numerical calculations is forced 
to pass through the first data point representing the average enthalpy of the 
gas entering the constant area constrictor of the plasma generator, and the 
axial distance is set equal to zero at this point. (A dashed line in parts 
(f) through (j) of the figures represents the analytical theory of refer- 
ence 1. In figures 10 through 14, the average enthalpy of the analytical 
theory was also made to match the average enthalpy of the gas entering the 
constant area constrictor in the plasma generator, as shown in part (h).) 

For the comparison with the experimental measurements within the 
0.635 -cm-diameter plasma generator (figs. 15-20), the numerical calculations 
were started in the stagnation region a short distance downstream from the tip 
of the cathode. It has been shown that a change in the starting conditions 
strongly affects the solutions only within a short distance downstream from 
the starting location. (See discussion under "Development of the Numerical 
Analysis.") Therefore, since the numerical solutions were started in the 
stagnation chamber, the solutions throughout the narrow constricted section 
were affected very little by the assumed starting conditions; so it was not 
necessary to match the average enthalpies at the inlet of the constant area 
constrictor when the numerical calculations were started in the stagnation 
chamber . 

The enthalpy at the start of the narrow constricted section was arbitrar- 
ily assumed to be zero for the analytical model of reference 1 in figures 15 
through 20. This theoretical model applies only to the constant area section, 
and for predictions from this theory when the inlet enthalpy is unknown, an 
arbitrary assumption of negligible enthalpy at the inlet has usually been 
employed. 

Values of gas radiance for the numerical calculations shown in figures 15 
through 20 were taken from reference 12 rather than reference 10 and are 
higher than the values given in reference 10. However, the radiation trans- 
port for these calculations is small compared to fhe thermal conduction; con- 
sequently, the solutions would be nearly the same if values of gas radiance 
from reference 10 had been used. 

The numerical calculations agree with the experimental measurements 
within a factor of 2 and illustrate the qualitative trends of the measure- 
ments. This agreement is within the agreement between the theoretical and 
the measured transport properties at high temperatures. 

The deviation between the numerical calculations and the experimental 
data for the 0 . 635 -cm constrictor is larger than for the 1 . 27 -cm constrictor. 
In the present numerical calculations the variation of the gas properties with 
pressure was extrapolated below 1 atmosphere. The large extrapolation to the 
low pressures for the comparison with the experimental results in the 0 . 635 -cm 
constrictor may be the cause of this larger discrepancy between the numerical 
solutions and the experimental data for the 0 . 635 -cm constrictor. 
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Arc -column characteristics were predicted within the accuracy with which 
the transport properties are known. The accuracy with which the theoretical 
model used in the numerical calculations represents the physical arc column 
cannot be determined from the above comparisons because of the uncertainties 
in the theoretical transport properties. 


STU US OF ARC CHARACTERISTICS WITH THE NUMERICAL PROGRAM 


The qualitative picture of the arc given by the numerical calculations 
(e.g., see fig. 10) shows that most of the gas is forced to the constrictor 
wall near the constrictor entrance; as the gas proceeds downstream, it is 
slowly ingested into the hot core of the arc (fig. 10(b)). The hottest region 
is near the constrictor entrance (fig. 10(a)). The energy flux densities 
(fig. 10(c)) .are highest at the constrictor exit rather than in the hot region 
near the entrance because near the entrance, most of the gas flows close to 
the cold wall rather than through the hot region. The radial velocity pro- 
files (fig. 10(d)) are approximately parabolic, similar to Poiseuille flow. 
(The velocity gradients are caused mainly by large density gradients rather 
than by viscous forces for this case. The effects of viscosity are illus- 
trated later in this section.) The distribution of momentum flux (fig. 10(e)) 
illustrates that near the constrictor inlet, momentum is convected to the 
constrictor walls causing the "ears" on the momentum flux profiles near the 
constrictor inlet. Farther downstream, as the gas is ingested into the arc 
core, the momentum is convected radially inward, causing high momentum flux 
in the center of the arc. (The ears downstream appear to be caused by large 
variations in gas viscosity since they disappear when the viscosity is set 
equal to zero.) 

A comparison of the analytical model of reference 1 with the numerical 
solutions indicates that many of the important qualitative trends to the con- 
stricted arc are predicted by the analytical model of reference 1 whenever 
the conduction heat losses are dominant. For example, figures ll(f)-(h) 
illustrate that the analytical theory gives a fair estimate of the enthalpy 
at the constrictor exit (i.e., at z/z Q = 0.1) and figures ll(i) and (j) 
illustrate that the theory gives the approximate values for the voltage gra- 
dient and wall heat transfer rate throughout the constrictor. 

For pressures over 1 atm within the 1. 27 -cm-diameter constrictor or for 
larger constrictors at atmospheric pressure, the radiation losses are dominant 
and the more complete numerical calculations are required for reasonable pre- 
dictions of the constricted arc behavior. Figures 10(h) and 10(j) illustrate 
that the theory of reference 1 predicts a higher mass average enthalpy and 
lower wall heat transfer rate than obtained from experimental measurements or 
from the numerical calculations for a case wherein radiation losses are large 
compared to conduction losses. 

The effects of various approximations in the theoretical model are 
illustrated in appendix A. A more complete discussion of these approximations 
and a comparison of the numerical solutions with several other theoretical 
models that contain different approximations are given in reference 7* 
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The following additional are characteristics were investigated with the 
numerical program: 

1. The effect of asymmetries at the inlet of the constrictor. 

2. The effect of radiation heat losses on the characteristics of the 
constricted thermal arc. 

3. The importance of viscous forces on the constricted thermal arc. 

4. The effect of radial convection within the constricted-arc column, 
including effects of transpiration -cooled constrictors. 

5. Behavior of the constricted arc at large dimensionless length. 

The numerical solutions obtained for these studies are illustrated in fig- 
ures 21 through 36 . Since the reader may wish to make many comparisons 
between the solutions in addition to those specifically mentioned herein, the 
solutions are presented completely and uniformly rather than for specific 
comparisons only. The solution of the simplified model in reference 1 is 
also included in parts (f) through (j) of these figures. 

Results of these studies are as follows: 

1. Any asymmetry that may be imposed on the axe at the constrictor 
inlet decays to a negligible value within a constrictor length equal to one- 
tenth the characteristic length, z G , as shown in figure 36 . 

2. The major effect of large radiation losses is to lower and flatten 
the enthalpy distribution (cf. figures 21-23 that are cases with negligible, 
moderate, and large ratios of radiation losses to conduction losses) . 

Note that the radiation does not appear to change the characteristic 
length greatly for the arc to approach an asymptotic state wherein the mass 
average enthalpy does not change with axial distance. The theory of refer- 
ence 1 predicts that the mass average enthalpy will approach 80 percent of 
the asymptotic value at z/z G = 0 . 1 . The solution in figure 23 indicates 
that the same may be true for arcs with high radiation. This result illus- 
trates that the spreading of the arc to the constrictor wall is due to thermal 
conduction; therefore, the theoretical characteristic length given in refer- 
ence 2 that has been modified to reflect a large dependence on radiation 
losses (modified length = z Q times fraction of heat losses due to thermal 
conduction) may underestimate the characteristic length considerably. From 
this modified characteristic length one would have predicted that the mass 
average enthalpy would have approached 80 percent of the final value at z/z D 
less than 0.03 for the case shown in figure 23 • 

3- The viscous forces are negligible compared to inertia forces in 
plasma generators of the size and with the gas flows normally used in wind 
tunnels (constrictor diameters over 1 cm and mass flows over l/2 g/s) . A 
comparison of figures 24 and 25 shows that neglecting the viscous forces in 
the numerical calculations does not change the solution appreciably. Further- 
more, when the viscous forces are negligible, the pressure drop through an 
aerodynamically choked constrictor is approximately half the inlet pressure. 



For smaller diameter constrictors at lower flow rates , viscous forces 
tend to flatten the velocity profiles. Figures 2 6 and 27 illustrate the 
effect of viscosity in a 0 . 635 -etn -diameter constrictor with a mass flow of 
only 0.227 g/s. 

4. The radial heat conduction is large compared to radial heat convec- 
tion within constricted arcs longer than 0 . 1 z Q (cf. figs. 28 and 29 that 
show the numerical calculations with and without the radial heat convection) . 
The solutions of the constricted arcs with part of the gas transpired through 
the constrictor wall do not differ appreciably from the solutions for the arcs 
where all of the gas enters at the constrictor inlet (provided the total gas 
flow rate is the same for both cases) , even at flow rates sufficiently large 
that radial convection is appreciable. Anderson (ref. 13 ) showed that radial 
convection can be made equal to radial conduction for gas flow rates greater 
than 2 Qn:k/cp times the constrictor length (or equivalently, with flow rates 
such that the constrictor length is less than 0.05 z Q ). Figures 30 through 33 
show that the enthalpy and energy flux at the constrictor exit are nearly the 
same in hydrogen or nitrogen arcs whether the gas is transpired through the 
walls or put in at the constrictor inlet. When all the gas enters at the con- 
strictor inlet, it rapidly spreads to the wall and then gradually returns 
toward the axis with a radial flow similar to the radial flow within the 
transpiration cooled constrictor. 

5- The enthalpy, mass flow, and energy flux distributions approach 
asymptotic values at large dimensionless lengths, but the velocity and momen- 
tum distributions continually increase in magnitude with length as shown in 
figure 3 ^* 


PREDICTIONS OF CONSTRICTED -ARC PLASMA. GENERATOR PERFORMANCE 


The numerical calculations indicate that the constricted arc may be used 
to produce gas flows with high energy flux density, with high enthalpy, and 
with high velocity. Figure 37 illustrates that an energy flux density of 
2600 kW/cm 2 may be obtained with hydrogen with a 0 . 635-cm-diameter constrictor 
that experiences a maximum wall heat transfer rate of only 2.4 kW/cm 2 , and 
figure 38 illustrates that an energy flux density of 1400 kW/cm 2 may be 
obtained with nitrogen in a 0 . 635 -cm-diameter constrictor that experiences a 
maximum wall heat transfer rate of 8 kW/cm.2. Figure 39 illustrates that with 
nitrogen, enthalpies over 5x10 s j/kg and velocities in excess of 9^00 m/s may 
be obtainable at the constrictor exit. (Note that the velocity at the exit is 
approximately sonic. A simple estimate of the exit velocity can be made from 
the enthalpy by noting that for nitrogen the sonic velocity is approximately 
equal to (0.l4 times enthalpy ) 1/2 for enthalpies over 2X10 7 j/kg as shown in 
fig. 40. ) The velocity at the exit of the divergent section of a plasma gen- 
erator such as shown in figure 1 should be more than twice as high as the 
velocity at the constrictor exit. Even if the gas were frozen and monatomic, 
the velocity at the nozzle exit would be nearly twice the velocity at the 
throat; and if the gas were not frozen or not monatomic, the velocity would 
be greater. Therefore, velocities in excess of 18,000 m/s may be obtainable 
at the nozzle exit . 


The energy flux density, enthalpy, and velocity given above cannot be 
obtained simultaneously; the conditions were chosen to optimize either the 
energy flux density or the enthalpy and velocity. (The methods to optimize 
the constrictor design to obtain the specific arc features are discussed in 
the next section.) Theoretical estimates of a more complete list of plasma 
properties that may be obtained simultaneously with a constricted-arc plasma 
generator similar to the one shown in figure 1 are given in table III for a 
single set of operating conditions. The theoretical estimates were based on 
a numerical solution for the constricted arc in the constant area portion of 
the constrictor; this numerical solution is shown in figure 4l. 

This table illustrates that the constricted arc should be capable of 
producing plasmas that interact appreciably with magnetic fields as required 
for experiments in magnetoplasmadynamics . 

Some interesting arc-column properties for a variety of operating 
conditions are shown in table II. 


DESIGN OPTIMIZATION OF CONSTRICTED -ARC PLASMA GENERATORS 


There is no constrictor size that is simultaneously optimum for all of 
the desirable performance characteristics of the constricted-arc plasma gen- 
erators, but the constrictor size, operating pressure, and current can be 
chosen to optimize any single performance characteristics. 

The features most often desired for these plasma generators are: 

(a) High enthalpy 

(b) High velocity 

(c) Broad uniform plasma stream for testing 

(d) High energy flux density 

(e) High electron density 

(f) High efficiency 

The limitations for obtaining the above features are the current that can be 
carried by the electrodes, the heat transfer rate that can be accommodated by 
the constrictor walls, and the voltage gradient that can be supported by the 
constrictor insulation. 

Any single characteristic above can be improved, and simultaneously some 
of the other performance characteristics reduced, if the constrictor size, 
the operating pressure, and the current are modified as follows: 

(a) High enthalpy 
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Maximum enthalpy can be achieved by operating with the maximum current 
that the electrodes can carry at a pressure sufficiently low that the radia- 
tion losses are negligible. The constrictor diameter should be just small 
enough that the heat transfer rate to the constrictor wall is the maximum that 
the wall can accommodate. (in the conduction dominant regime, the enthalpy is 
approximately proportional to i/R, as shown in ref. 1.) Furthermore, the 
maximum local enthalpy occurs on the center line near the cathode, so for 
maximum enthalpy, the constrictor length should be short. 

(b) High velocity 

The highest velocities at the constrictor exit can be obtained with 
aerodynamically choked constrictors that produce the highest enthalpies; so 
the optimization given above for high enthalpy also applies to high velocity. 

(c) Broad, uniform plasma stream for testing 

The uniform portion of the test stream can be increased by increasing 
the pressure until the radiation losses flatten the profile (cf. figs. 4l(a) 
and 42(a)) and by lengthening the constrictor until the enthalpy profile 
becomes asymptotic. The diameter of the constrictor should be just large 
enough that the heat transfer rate to the wall is the maximum that the con- 
strictor wall can accommodate, (in the regime wherein radiation losses are 
dominant the wall heat transfer rate increases with increasing diameter.) 

The current should be sufficiently great to heat the gas to a temperature that 
will produce high radiation; that is, to approximately 15,000 K (see 
figs. 2(q) and 2(w)). Too high a current heats the gas in the center of the 
arc to a temperature sufficiently high that the radiation losses decrease 
with increasing temperature, and, despite high pressures, causes peaking in 
the enthalpy profile as shown in figure 43* 

(d) High energy flux density 

The energy flux density can be increased by decreasing the constrictor 
diameter, increasing the constrictor length until an asymptotic profile is 
obtained, and increasing the pressure. The current should be just high 
enough that the heat transfer rate to the constrictor wall is the maximum that 
the wall can accommodate. Higher energy fluxes can be obtained with the 
lighter gases because they radiate the least. (Compare the two cases shown 
in figs. 37 and 38 wherein the energy flux density with hydrogen is approxi- 
mately twice that with nitrogen, even though the heat transfer rate to the 
constrictor wall for the nitrogen is greater than that with the hydrogen.) 

(e) High electron densities 

The number density of electrons can be increased by decreasing the 
diameter and increasing the length until the profile becomes asymptotic. The 
current should be just high enough that the enthalpy of the gas is that for 
which the maximum electron density occurs at constant pressure, and the pres- 
sure should be just high enough that the wall heat transfer rate is the 
maximum that the wall can accommodate. 
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(f) Efficiency 

The efficiency can be increased by operating at low pressure with short 
constrictors. At low pressure, the efficiency is not strongly influenced by 
either the constrictor diameter or by the current. If high pressure operation 
is mandatory, then for some cases, contrary to all of the simplified theories, 
the efficiency may be improved while maintaining the same mass average 
enthalpy at the exit by lengthening the constrictor and reducing the current. 
(Compare figs. 42 and 44 that illustrate that the efficiency was increased 
from 35 to 37 percent by increasing the length and decreasing the current 
while maintaining the same mass average enthalpy at the exit.) 


CONCLUSIONS 


The equations that govern the energy and mass transport in the con- 
stricted arc have been solved numerically. The numerical program uses real 
equilibrium gas properties and solves simultaneously the energy, momentum, 
and continuity equations. Axial conduction, radial pressure gradients, and 
radial voltage gradients are neglected. From given initial and boundary con- 
ditions the numerical program solves for the local state properties and 
velocity of the gas within the constricted thermal arc. From these properties 
all other characteristics of interest (i.e., voltage gradient and wall heat 
transfer rates) can be obtained. These solutions give sufficient detail to 
evaluate other approximate solutions, and the computing time (approximately 
2 min) is sufficiently short that the programs may be used directly to obtain 
design criteria for plasma generators. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif., 94035, Oct. 13, 1966 
129-01-02-01-00 -21 
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APPENDIX A 


NUMERICAL SOLUTIONS FOR EVALUATING APPROXIMATIONS OF THE 
SIMPLIFIED THEORETICAL MODELS OF THE CONSTRICTED ARC 
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Figure Al.- The complete numerical solution without the simplifying approximations. 
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Figure A2.- Numerical solution with the approximation that mass flux density is constant. 
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Figure A^. - Numerical solution with the idealized gas of reference 1. 
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(f) Center -line enthalpy. 



k 


(g) Space-average enthalpy. 
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(c) Energy flux density. 


Figure A5-- Numerical solution with the assumption that the initial enthalpy distribution is a 

Bessel function. 
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Figure A 6.- The simplified theoretical model of reference 1. (This model contains all of the 

approximations in figures A2 through A5- ) 
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Figure AT.- The model of reference 1 without the approximation that the mass-flow density is constant. 







(d) Center -line enthalpy. 



(e) Space -average enthalpy. 



.02 .04 .06 .08 .10 

z/z 0 
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Figure A9. - The model of reference 1 without the assumption that the initial enthalpy distribution is 

a Bessel function. 




(e) Space-average enthalpy. 



J 1 1 1 1 I 

.02 .04 .06 .08 .10 .12 

z/z 0 

(f) Mass -average enthalpy. 


Figure A9- 






38 



APPENDIX B 


FORTRAN PROGRAM FOR THE SOLUTION OF THE AXISYMMETRIC CONSTRICTED ARC 

WITH AN AXIAL FLOW OF GAS 


This appendix contains the Fortran II programs for the numerical solution 
of the axisymmetric constricted arc with an axial flow of gas. The data 
required for input into these programs and a legend for the Fortran variables 
is also included. 

The function of each subroutine is described below. 


BOUNDC provides for the specification of the boundary conditions for the 

main program 

STATEP evaluates the state properties (except density) for all of the 
mesh points at each axial station 

WDOT evaluates the density for all of the mesh points at each axial 

station and calculates the mass flow rate by integration of the 
mass flux over the constrictor cross-sectional area 


ITER provides the iteration of pressure to obtain the correct mass flow 

in the subsonic portion of the nozzle 


ITERS provides the iteration of pressure to obtain the correct mass flow 

in the supersonic portion of the nozzle 

MOM calculates the velocity for all mesh points at the next axial 

station from the momentum equation 


ENERGY calculates the enthalpy for all mesh points at the next axial 

station from the energy equation 

OUTPT prints and writes on magnetic tape the results of the main program 
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INPUT DATA REQUIRED 


Card number 

Data or Fortran variable 

Format 

1 

NFILE 

14 

2 

NTAPE 

14 

3 

KMAX 

14 

4 

KINC 

14 

5 

AMPS 

ElO .3 

6 

WS 

EIO .3 

7 

TRCL 

E 10.3 

8 

P(l) 

ElO .3 

9 

DIAM(l) 

E 10.3 

10 

THETA 

ElO .3 

11 

HWALL(l) 

ElO .3 

12 

NMESH 

14 

13 

FZO 

E 10.3 

14 

EX 

E 10.3 

15 

EXX 

ElO .3 

16 

EPS 

E 10.3 

17 through 17 + N 

H(l) through H(NMESH) 

5 EIO .3 

l8 + N through l8 + 2N 

U(l) through U(NMESH) 

5E10.3 


where 

N = MOD[ (NMESH - l)/5l 

These Fortran variables are described in the following legend under the 
variables common to all of the programs or under the variables for program 
number HT0702. 


VARIABLES IKED IN THE PROGRAMS FOR THE NUMERICAL SOLUTIONS OF THE 
SYMMETRIC CONSTRICTED THERMAL ARC WITH AN AXIAL FLOW OF GAS 

Variables Common to All Programs 

Variable Name Description 

AMPS total current carried by the constrictor, A 

DIAM diameter of the constrictor (an array) , m 

DIA diameter of the constrictor, m 

DP pressure drop between axial stations, N/m s 

DRDR square of the incremental radial distance, DR, m 2 

DR incremental radial distance, m 

DW discrepancy between the mass flow rate at this axial 

station and the initial mass flow rate, kg/s 


Variable 

DZ 

EPS 

E 

EX 

EXX 

FZO 

HAVE 

HRAVE 

H 

HWALL 

KINC 

KMAX 

K 

LOC 

L 

M 

NCHOKE 

NERR 

NFILE 

m 


name _ Descripti on 

incremental axial distance, m 

maximum allowable relative discrepancy of the mass flow 
rate 

axial voltage gradient, V/m 

factor by which the axial incremental distance is increased 
for the next axial station 

stability factor that limits maximum size of DZ in order 
to maintain stability 

length of the first axial increment divided by the 
characteristic arc length, z Q 

space average of the enthalpy, — J h dA, j/kg 

A d a 

average energy density, “ J^ph dA, j/m 3 
local enthalpy, j/kg 

enthalpy of the gas at the constrictor wall, j/kg 
number of axial stations between printout of results 
maximum number of axial stations to be calculated 
axial station number 

number of iterations required to satisfy the continuity 
equation (i.e., the number of iterations required such 
that the flow rate at this axial station is suff iciently 
near the initial flow rate) 

relative axial station number 

second relative axial station number 

flag indicating that the flow is choked 

error flag for the gas property subroutines 

file number on the magnetic tape used to store the 
solutions 

extra variable not used 
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Variable name 
NMESH 

NNN 

NTAPE 

PHI 

PHIW 

P 

QP 

Q 

ECAP 

RHOAV 

RHO 

RHOU 

REOU 

R 

RUH 

SIGMA 

THETA 

TRCL 

U 

vise 

¥ 


Description 

number of radial increments from the center line to the 
constrictor wall 

flag indicating this is the first iteration on pressure 
drop (zero indicates first iteration, 1 indicates all 
iterations thereafter) 

magnetic tape on which the solution is stored 

local thermal conductivity potential, W/m 

thermal conductivity potential of the gas at the constric 
tor wall, W/m 

static pressure at this axial station, N/m 2 

heat transfer rate to the constrictor wall that is due to 
radiation, W/m 2 

heat transfer rate to the constrictor wall that is due to 
thermal conduction, W/m 2 

local radiation, W/m 3 

space average density, kg/m 3 

local density, kg/m 3 

local product of density and velocity at this station, 
kg/sm 2 

local product of density times velocity at the previous 
axial station, kg/sm 2 

local radius, m 

total energy flux at this axial station, W/m 2 
local electrical conductivity, l/R-m 

half-angle of divergence of the supersonic nozzle, deg 

mass flux density of constrictor transpiration cooling, 
kg/sm 2 

local velocity, m/s 
local viscosity, Ns/m 2 

mass flow rate at the initial axial station, kg/s 


b2 



Variable name 


Description 

WW mass flow rate at this axial station, kg/s 

Z axial distance of this axial station from the initial 

axial station, m 


DZMAX 

FNMESH 

KC 

KI 

KMESHP 

m 

ms 

PRAD 

PRES 

Q3? 

RUHA. 

THETAR 

ZC 


Variables Local to Program HT-O7OI 

upper limit of the axial increment distance (the incre- 
mental axial distance is never allowed to exceed this 
value in order to keep the solution stable) , m 

number of radial increments from the center line to the 
constrictor wall (floating number) 

axial station number at which aerodynamic choking occurs 

index for the main loop within the program 

radial index for the constrictor wall 

index to indicate the axial stations for which the results 
should be printed and stored 

flag to indicate that the flow is supersonic 

percentage of the heat flux to the constrictor wall at 
this axial station that is due to radiation 

local static pressure, atm 

total heat transfer rate to the constrictor wall at this 
axial station, W/m 2 

mass average energy at this axial station 
j/kg 

half-angle of divergence of the supersonic nozzle, radians 

axial distance from station 1 for which the flow chokes, m 


1 r» 

, puh dA, 

9 m Ja 


CWF 

RHOA 


Variables for Program HT-O7O2 

scaling factor that changes the magnitude of the velocity 
profile to obtain the desired initial flow rate 

space average of the density at this axial station, kg/m 3 
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Variable name 
THEPAR 

WS 

ZO 


Description 


half-angle of divergence of the supersonic nozzle, radians 
mass flow rate through the constrictor 
characteristic arc length, z Q , m 


DA 

DHA. 

DHRA 

DQR 

DBS 

FJ 

HA. 

HRA 

J 

HMESHP 

PRES 

RDR 

RW 

SS 

SW 

VW 


Variables for Program HT-0703 
incremental cross-sectional area, m 2 

product of the local enthalpy and the incremental area, 
Jm 2 /kg 

product of the local enthalpy, the local density, and the 
incremental area, j/m 

product of the local radiation and the incremental area, 
W/m 

product of the local electrical conductivity and the 
incremental area, m/fi 

floating index for the radial position 

running total of the product of the enthalpy and the area, 
Jm 2 /kg 

running total of the product of enthalpy, density, and the 
area, j/m 

index for the radial position 
radial index for the constrictor wall 
local static pressure, atm 

product of the local radius and incremental radial 
distance, m 2 

dummy variable not used 

running total of the product of electrical conductivity 
and area, m/n 

electrical conductivity at the constrictor wall, l/n-m 
viscosity at the constrictor wall, Ns/m 2 
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Variables Local to Program HT -0704 


Variable 

DA 

DBA 

DRUH 

J 

PRES 

HA 

RDR 

DDP 

IND 

N 

AH 

DA 

DRUTP 

DR 17 T 

GA 

GB 


name Description 

incremental cross-sectional area, m 2 

product of the incremental cross-sectional area and the 
local mass flux, kg/s 

product of the incremental cross-sectional area and the 
local energy flux, W 

index for the radial incremental distance 
local static pressure, atm 
running sum of the mass flux, kg/s 

product of the radius and the incremental radial distance, 
m 2 


Variables Local to Programs HT-O705 and 0706 

amount that the pressure drop is changed for each itera- 
tion at this axial station, N/m 2 

flag used in adjusting the pressure drop during the itera- 
tions at this station 

index for the iterations 


Variables Local to Program HT-O7O7 

cross-sectional area of the previous station divided by 
the cross-sectional area of this station 

incremental cross-sectional area, m 2 

radial mass flux from this volume increment toward the 
center of the column, kg/sm 2 

radial mass flux toward the center of the column into this 
volume increment, kg/sm 2 

velocity of the gas associated with the radial mass flux, 
DRUTP, m/s 

velocity of the gas associated with the radial mass flux, 
DROT, m/s 
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Variable 

NMESHP 

RADCON 

RCP1 

RCP2 

RCP3 

RCP4 

RDR 

RUPREV 

VL 

AR 

CL 

DA 

DRUTP 

DRUT 

FJ 

GA 

GB 

J 

46 


name 


_ Descript ion 

index for the radial incremental distance that designates 
the constrictor wall 

momentum convected radially into this volume increment , 
kg/ms 2 

momentum convected radially into the volume element at the 
previous axial station, kg/ms 2 

momentum convected radially into the volume increment two 
axial stations ahead of this one, kg/ms 2 

momentum convected radially into the volume increment 
three axial stations ahead of this one, kg/ms 2 

momentum convected radially into the volume increment four 
axial stations ahead of this one, kg/ms 2 

product of the radius and the radial incremental distance, 
m 2 

local mass flux at the previous axial station, kg/sm 2 

viscous losses for this volume increment, N/m 2 


Variables Local to Program HT-0708 

cross-sectional area of the previous station divided by 
the cross-sectional area of this station 

conduction losses from this incremental volume, W/m 3 

incremental cross-sectional area, m 2 

radial mass flux from this volume increment toward the 
center of the column, kg/sm 2 

radial mass flux toward the center of the column into this 
volume increment, kg/sm 2 

index for the radial incremental distance (floating) 

energy of the gas associated with the radial mass flow, 
DRUT , j/kg 

energy of the gas associated with the radial mass flow, 
DRUTP, J/kg 

index for the radial incremental distance 
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Variable 

OH 

RADCON 

RCP1 

RCP2 

RCP3 

RCP4 

RDR 

RL 

TE 

TEP 

PRAD 

PRES 

QT 

RUHA 


name 


_ Description 

ohmic heating within this volume element. 


W/m a 


energy convected radially into this volume increment , 
kg/m s 2 

energy convected radially into the volume element at the 
previous axial station, kg/ms 2 

energy convected radially into the volume increment two 
axial stations previous to this one, kg/ms 2 


energy convected radially into the volume increment three 
axial stations previous to this one, kg/ms 2 

energy convected radially into the volume increment four 
axial stations previous to this one, kg/ms 2 


product of the radius and the radial incremental distance, 
m 

local radiation, W/m 3 

local kinetic energy (an array), j/kg 
local kinetic energy, j/kg 


Variables Local to Program HT-O 7 O 9 

percentage of the constrictor wall heat flux that is due 
to radiation 

local static pressure, atm 

total heat transfer rate to the constrictor wall, W/m 2 

mass average energy at this axial station, ^ J\puh 
j/kg m A 





n n 


CHT0701 VAL WATSON DC ARC - AXIAL FLOW - REAL GAS - 2D 

COMMON AMPS * DIAM* DIA* DP* DRDR * DR* DW* DZ * EPS* E 
1 * EX* EXX* FZO * HAVE* HRAVE* H* HWALL* K* KINC* KMAX » LOC* L 

1 ♦ M» NCHOKE * NERR. NFlLE* NK* NMESH • NNN » NTAPE* PHI* PHIW 

1 * P* QR* Q* RCAP * RHOAV* RHO* RHOU* RROU* R* RUH 

1 * SIGMA* THETA* TRCL * U» VISC* W* WW* Z 

DIMENSION D I AM ( 2000 ) * P(2000)* E ( 2 0 0 0 ,) • HWALL(2000) 

1 * PHI f 100 ) * SIGMA ( 100 ) * RCAP ( 100 ) ♦ VISC(IOO)* RHO(IOO) 

I * R ( 100 ) * RHOUUOO)* RROU ( 100 ) 

1 * H ( 2 * 100 ) * U ( 2*100 ) 

CALL CRISISCZZ ♦ AA * AMPS *Z ) 

98 CONTINUE 
NSS » 0 
C 

C SET INITIAL CONDITIONS AND COMPUTE FIRST AXIAL STEP 
CALL BOUNDC 
FNMESH « NMESH 
C 

REWIND 8 

C MAIN LOOP FOR COMPUTING EACH AXIAL STEP 
DO 6 K I s 3 *KMAX 
K * K +1 
NN = NN - 1 
C 

C MAINTAIN AXIAL STEP SIZE LESS than STEP SIZE FOR INSTABILITY 

DZMAX = ( (DIA#DIA/4.0)*RH0U(2)*H(M*2) /(PHI I 2 ) *FNMESH*FNMESH)) *EXX 
IF(DZ-DZMAX) 40*42*42 
40 DZ » EX*DZ 

42 CONTINUE 

Z = Z + DZ 

START DIVERGENCE OF NOZZLE AFTER CHOKING 
I F ( NSS ) 58*58*52 
52 IFtZ-ZC) 58*56*56 

56 DIAM(K) = DIAM(KC) +2. 0* ( Z-ZC ) *TANF ( THETAR ) 

58 CONTINUE 

DIA = DIAM(K) 

DR » DIA/(2.0#FLOATF(NMESH-m 
DRDR a DR*DR 
C 

C INCREASE IN FLOW RATE FROM TRANSPIRATION COOLING 
W a w + DZ*3.1416*DIA#TRCL 
C 

C ALTERNATING STORAGE LOCATION FOR AXIAL STATIONS 
M = 1 
L = L + 1 
IFI3-L) 1*1*2 



no on 


i 


M = 2 
L = 1 

2 CONTINUE 

EVALUATION OF THE ENTHALPY AT NEXT AXIAL STATION FROM ENERGY EQUATION 
CALL ENERGY 

CHECK FOR SUPERSONIC OR SUBSONIC FLOW 
IF(NSS) 60*60*62 

C CALCULATION OF VELOCITY AT NEXT STA THRU ITERATION - SUBSONIC FLOW 
60 CALL ITER 
GO TO 68 

62 IF(Z-ZC) 64*66*66 
64 CALL ITER 
GO TO 68 
C 

C CALCULATION OF VELOCITY AT NEXT STA THRU ITERATION - SUPERSONIC FLOW 
66 CALL ITERS 
68 CONTINUE 
C 

C CHECK FOR CHOKED FLOW 

C IF CHOKED AND SUBSONIC* START DIVERGING NOZZLE 
C IF CHOKED AND SUPERSONIC* GIVE ERROR READING AND EXIT 
I F ( NCHOKE ) 4*4*70 
70 IF(NSS) 72*72*3 

72 NSS = 1 

KC = K - 1 
ZC = Z - DZ/2.0 

CALL SKIP (“1 *NTAPE) 

READ TAPE NTAPE* K, DIAM(K)* Z> AMPS* E ( K ) * W»QT *HWALL (K) 

1 *PRAD* NMESH * PRES * FZO *LOC ♦ EPS »DW* DP* DZ 

2 * HAVE* RUHA*HRAVE 
P ( K ) = PRES*1.013E5 

READ TAPE NTAPE* (R(J1* H(M*J)* U(M*J>* RHOU(J)* J=1*NMESH) 

NMESHP = NMESH + 1 
U ( M ,NMESHP ) = 0.0 - U ( M * NMESH ) 

DO 73 J=1«NMESH 

73 RROUf J)=RH0U( J) 

NN=KC-K 

CALL SKIPIl. NTAPE) 

CALL STATEP 
CALL WDOT 
GO TO 6 

3 WRITE OUTPUT TAPE 6* 202* K* DW* U(M*2) 

GO TO 8 

4 CONTINUE 


<o 



o o o o 


CALL STATEP 




VJl 


IF PRESSURE TOO LOW FOR GAS TABLES* EXIT 
IF ( P ( K ) - 0.1E5) 74*74*76 
74 CALL OUTPT 

GO TO 8 
76 CONTINUE 

WRITE OUT VALUES FOR EVERY (KINC)TH AXIAL STATION 
IF(NN) 5*5*6 

5 NN = KINC 
CALL OUTPT 

6 CONTINUE 
8 CONTINUE 

REWIND NTAPE 

202 FORMAT { 1H0* 18HFLOW CHOKED AT K = * 14* 10X* 

1 20HFLOW RATE ERROR IS ♦ 2PF12»7* 10H PERCENT » 

2 14 HCL VELOCITY = *0PF10.1* 8H M/SEC ) 

GO TO 98 

END 




\ 
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CHT0702 VAL WATSON SUBROUTINE BOUNDC FOR HT0720 
SUBROUTINE BOUNDC 

COMMON AMPS » DIAM* DIA* DP* DRDR* DR* DW* DZ* EPS* E 
1 * EX* EXX* FZO* HAVE* HRAVE H* HWALL* K* KINC* KMAX * 

1 * M* NCHOKE • NERR* NFILE* NX* NMESH * NNN * NTAPE • PHI* 

1 ♦ P* QR* Q* RCAP • RHOAV* RHO* RHOU* RROU* R* RUH 

1 » SIGMA. THETA* TRCL ♦ U-» VISC* W* WW* Z 

DIMENSION D I AM ( 2000 ) » P ( 2000 ) * E(2000)» HWALL(2000) 

1 * PHI { 100 ) » S I GMA ( 100 ) * RCAP ( 100 ) * V I SC ( 100 ) * RHO ( 

1 * R(IOO). RHOU ( 100) * RROU ( 100 ) 

1 * H ( 2 * 1 00 ) * U ( 2 * 100 ) 

CALL CRISIS (ZZ.AA ) 


LOC* L 
PHIW 


100 ) 


SET UP MAGNETIC TAPES 

READ INPUT TAPE 5* 100* NFILE * NTAPE 
CALL LOCATE (NFILE*NTAPE) 

SET MAX ALLOWABLE AXIAL STATIONS AND INTERVAL BETWEEN PRINTOUT 
READ INPUT TAPE 5* 100* KMAX* KINC 

SET INITIAL VALUES OF CURRENT* FLOW RATES* AND PRESSURE 
READ INPUT TAPE 5* 101* AMPS* WS* TRCL* P ( 1 ) 

P ( 1 ) = P ( 1 ) *1*013E5 

SET UP THE DIA AS A FUNCTION OF AXIAL DIST 
DIAMETER IS CONSTANT IN THIS CASE 

itheta is the half angle of divergence after choking occurs) 

READ INPUT TAPE 5. 101* DIAM(l) 

DO 300 ID * 2*KMAX 
300 D I AM ( ID) = DIAM(l) 

READ INPUT TAPE 5* 101* THETA 
THETAR = THETA* 2* 0*3* 1416/360*0 

SET THE WALL TEMP EQUAL TO A FUNCTION OF AXIAL DISTANCE 
WALL TEMP IS CONSTANT IN THIS CASE 
READ INPUT TAPE 5* 101* HWALL ( 1 ) 

DO 400 IW = 2 *KMAX 

HWALL (IW) = HWALL ( 1 ) 

THE RADIAL MESH SIZE* THE RELATIVE AXIAL INCREMENT SIZES* 

THE STABILITY FACTOR ( RATIO OF MESH SIZES)* AND THE FLOW RATE ACCURACY 
READ INPUT TAPE 5*100* NMESH 
READ INPUT TAPE 5* 101* FZO* EX* EXX* EPS 

THE INITIAL ENTHALPY AND VELOCITY RADIAL PROFILES AT FIRST STATION 
READ INPUT TAPE 5*102* (H(1*J)» J=1*NMESH) 

READ INPUT TAPE 5.102* (Ufl.J)* J=1*NMESH) 


400 
SET 

C 

C SET 



n n 


r 


c 

C EVALUATE THE REMAINING PROPERTIES AT THE FIRST AXIAL STATION 
K = 1 

DIA = DIAM(K) 

DR = D I A / ( 2 »0*FLOATF ( NMESH-1 ) ) 

DRDR = DR#DR 
L = 1 
M = 1 
Z = 0.0 
CALL STATEP 
CALL WDOT 

ADJUSTMENT FOR PROPER FLOW RATE 
CWF = WS/WW 
DO 2 J=1 »NMESH 

U ( 1 * J ) = CWF*U ( 1 * J ) 

2 U ( 2 * J ) = U(l.J) 

CALL WDOT 
W h WW 
C 

C SET THE INITIAL AXIAL INCREMENTAL DISTANCE EQUAL TO FZO*CHARACT. LENGTH 
ZO = W*H(1*2)/(PHI(2)*3.1416) 

DZ = FZO*ZO 
CALL OUTPT 
C 

C CALCULATE THE PROPERTIES FOR THE SECOND AXIAL STATION 
M = 2 
K = 2 

DIA = DIAM(K) 

DR = DIA/(2.0*FLOATF(NMESH-1M 
DRDR = DR*DR 
DP = 0.0 

P(K) = P(K-l) + DP 
CALL ENERGY 
RHOA = RHOAV 
CALL WDOT 

DP =W#W* ( RHOAV— RHOA ) / ( ( (RHOA*3.1416*DIA*DIA) /4»0 )#*2 ) 

CALL ITER 
Z = Z + DZ 
CALL STATEP 
CALL OUTPT 

100 FORMAT ( 14 ) 

101 FORMAT ( E10. 3 ) 

102 FORMAT ( 5E10 .3 ) 

RETURN 

END 



Ui 

U) 
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CHT0703 VAL WATSON SUBROUTINE STATEP FOR HT0720 
SUBROUTINE STATEP 

COMMON AMPS * DIAM* DI A* DP* DRDR * DR* DW* DZ» 

1 * EX* EXX* FZO * HAVE* HRAVE* H. HWALL* K. * 

1 * M » NCHOKE * NERR » NFILE* NK* NMESH * NNN » 

1 * P* QR* Q* RCAP ♦ RHOAV* RHO* RHOU* RROU* 

1 


1 

EPS* E 

KINC* KMAX * LOC* L 
NTAPE* PHI* PH I W 
R* RUH 


* SIGMA* THETA* TRCL* U* VISC* W* WW» Z 
DIMENSION D I AM( 2000 ) * PI2000)* E ( 2000 ) * HWALL { 2000 ) 

1 * PHI ( 100 ) * S I GMA (100) * RCAP ( 100 ) * VISC(IOO)* RHOIIOO 

1 * R ( 100 ) * RHOU (100) • RROU ( 100 ) 

1 * H ( 2 * 100 ) * U ( 2*100 ) 

CALL CRISIS ( ZZ * AA ) 


EVALUATION OF THE GAS PROPERTIES AT THE WALL TEMPERATURE 
PRES = P(K) /1.013E5 

CALL NTABtPRES * HWALL ( K ) * PHIW* SW* RW« VW* 8* NERR) 
NERR = NERR 

HA = 0.0 

SS = 0.0 

HRA * 0.0 

QR = 0.0 

DO 30 J=1 *NMESH 


EVALUATION OF THE GAS PROPERTIES AT EACH RADIAL STATION 

CALL NTAB(PRES*H(M* J) *PHI ( J ) *SIGMA( J) *RCAP( J) * VI SC ( J ) * 8* NERR) 
NERR = NERR 
IF (NERR) 10*10*40 
10 CONTINUE 

IF (J-l) 20*30*20 
20 CONTINUE 

FJ = J 

R ( J ) = ( FJ-1 . 5 ) *DR 
RDR = R(J)*DR 
DA = 6 • 2832#RDR 
DHA = DA*H(M*J) 

DSS = DA*SIGMA ( J) 

DHRA = DHA*RHO(J) 

DQR = DA*RCAP(J) 

HA = HA + DHA 
SS = SS + DSS 
HRA = HRA + DHRA 
QR = QR + DQR 
30 CONTINUE 

NMESHP = NMESH + 1 

PHI ( NMESHP ) = 2.0*PHIW - PHI (NMESH) 

R ( NMESHP ) = DIA/2.0 
H(M .NMESHP) = HWALL (K) 



r 

I 


VISC ( NMESHP ) = 2 • 0*VW - VISC(NMESH) 

C 

C CALCULATION OF THE VOLTAGE GRADIENT* AVE ENTHALPY * AND HEAT FLUXES 
E ( K ) = A M.PS/SS 

HAVE = HA/(3il416*DIA*DIA/4.0) 

HRAVE= HRA/( 3.1416*DIA*DIA/4.0) 

QR = QR / ( 3 • 1416*D 1 A ) 

Q = 2« 0* ( PHI ( NMESH ) - PHIW)/DR 
40 CONTINUE 
RETURN 
END 



Ul 

VJ1 



vn CHT0704 VAL WATSON SUBROUTINE WDOT FOR HT0720 

^ SUBROUTINE WDOT 

COMMON AMPS * DIAM* DIA* DP* DRDR ♦ DR* DW * DZ* EPS* E 
1 * EX* EXX* FZO * HAVE* HRAVE* H* HWALL* K* KINC* KMAX * LOC* 

1 . M* NCHOKE » NERR* NFILE* NK. NMESH * NNN * NTAPE • PHI* PHIW 

1 * P* QR* Q* RCAP * RHOAV* RHO* RHOU* RROU* R* RUH 

1 * SIGMA* THETA* TRCL ♦ U* VISC* W* WW* Z 

DIMENSION D I AM ( 2000 ) • P ( 200 0 J • E(2000)* HWALL ( 2000 ) 

1 » PHI (100) * S I GMA (100)* RCAP ( 100 ) * VISC(IOO)* RHO(IOO) 

1 * R ( 100 ) • RHOU (100) ♦ RROU (100) 

1 * H ( 2 * 1 00 ) » U ( 2 » 100 ) 

CALL CRISIS (ZZ.AA ) 

C 

C EVALUATION OF THE FLOW RATE* AVERAGE DENSITY* AND ENERGY FLUX 
WW = 0.0 
RA = 0.0 
RUH = 0.0 

PRES = PHO/1.013E5 
DO 30 J=1 *NMESH 

CALL NRHO ( PRES * H(M*JJ* RHO(J>* 8* NERR) 

NERR = NERR 

RHOU(J) = RHO(JI*U(M* J) 

IF(J-l) 20*30*20 
20 CONTINUE 

RDR = R ( J ) *DR 
DA = 6 » 2832*RDR 
WW = WW + RHOU ( J) #DA 
DRA = DA*RHO(J) 

RA = RA + DRA 

DRUH = DA*H(M*J)*RHOU( J) 

RUH = RUH + DRUH 

30 CONTINUE 

RHOAV = RA/(3.1A16*DIA*DIA/4.0) 

RETURN 

END 




nnon 
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CHT0705 VAL WATSON SUBROUTINE ITER FOR HT0720 
SUBROUTINE ITER 

COMMON AMPS* D I AM * DIA. DP* DRDR * DR* DW* DZ * EPS* E 
1 * EX* EXX* FZO* HAVE* HRAVE* H* HWALL* K* KIN C* KMAX » LOC* L 

1 * M* NCHOKE* NERR* NFlLE* NK* NMESH * NNN • NTAPE » PHI* PHIW 

1 » P* OR* Q* RCAP * RHOAV * RHO* RHOU* RROU* R* RUH 

1 * SIGMA* THETA* TRCL* U* VISC* W» WW* Z 

DIMENSION D I AM ( 2000 ) * P f 2000 ) * E ( 2000 ) • HWALL(2000) 

1 * PHI ( 100 ) » SIGMA(IOO) » RCAP ( 100 ) * VISCtlOO)* RHO(IOO) 

1 » R ( 100 ) » RHOU (100) * RROUt 100 ) 

1 * H ( 2 * 100 ) ♦ U ( 2 * 100 ) 

CALL CRISIS (ZZ * AA ) 

ITERATION TO CALCULATE THE VELOCITY AT THE NEXT AXIAL STATION - SUBSONIC 
VELOCITY IS FROM MOMENTUM EQUATION 
ITERATE UNTIL THE MASS FLOW IS CONSERVED 
IND = 0 
NCHOKE = 0 
DDP = DP 

P ( K ) = PtK-1) + DP 
LOC = 0 
NNN = 0 
DO 15 N — 1 *50 
LOC * N 
CALL MOM 
NNN = 1 
CALL WDOT 
DW * (WW - W)/W 
IF( ABSF ( DW 1 - EPS) 20*1*1 
1 IF(DW) 3*20*9 

3 IF(IND) 7*7*5 

5 DDP = DDP/2.0 

7 DP = DP + DDP 

P ( K) = P(K-l) + DP 
GO TO 15 

9 DDP = DDP/2*0 

DP = DP - DDP 
P ( K ) = PIK-1 ) + DP 
IND = 1 

15 CONTINUE 
NCHOKE = 1 
20 CONTINUE 
RETURN 
END 



Ui 



noon 


CHT0706 VAL WATSON SUBROUTINE ITERS FOR HT0720 
SUBROUTINE ITERS 

COMMON AMPS » D I AM » DIA* DP* DRDR * DR* DW* DZ * EPS* E 
1 » EX* EXX* FZO * HAVE* HRAVE* H* HWALL* K • KINC* KMAX * LOC* L 

1 * M* NCHOKE * NERR* NFILE* NK* NMESH • NNN * NTAPE * PHI* PHIW 

1 » P* QR » Q* RCAP • RHOAV * RHO* RHOU* RROU* R* RUH 

1 * SIGMA* THETA* TRCL * U* VISC* W* WW* Z 

DIMENSION D I AM ( 2000 ) • P(2000)* E(2000), HWALL ( 2000 ) 

1 * PHI ( 100) * S IGMA ( 100 ) * RCAP ( 100 ) * VISC(IOO)* RHO(IOO) 

1 * R ( 100 ) * RHOU (100) * RROU ( 100 ) 

1 * H ( 2 * 1 00 ) * U ( 2 * 1 00 ) 

CALL CRI SIS (ZZ *AAA ) 



ITERATION TO CALCULATE THE VELOCITY AT THE NEXT AXIAL STATION - SUPERSONIC 
VELOCITY IS FROM MOMENTUM EQUATION 
ITERATE UNTIL THE MASS FLOW IS CONSERVED 
IND = 0 
NCHOKE = 0 
DDP = DP 

P ( K ) = P { K— 1 ) + DP 
LOC = 0 
NNN = 0 
DO 15 N=1 *30 
LOC = N 
CALL MOM 
NNN = 1 
CALL WDOT 
DW = (WW - W)/W 
IFlABSF(DW) - EPS) 20*1*1 
1 IF(DW) 9*20*3 

3 I F ( IND) 7*7*5 

5 DDP = DDP/2.0 

7 DP = DP + DDP 

P ( K ) = P(K-l) + DP 
GO TO 15 

9 DDP o DDP/2.0 

DP = DP - DDP 
P(K) = P(K-l) + DP 
IND = 1 

15 CONTINUE 

NCHOKE = 1 
20 CONTINUE 
RETURN 
END 
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CHT0707 VAL WATSON SUBROUTINE MOM FOR HT0720 
SUBROUTINE MOM 

COMMON AMPS » DIAM* DIA* DP* DRDR * DR* DW* DZ* EPS* E 
1 * EX* EXX. FZO » HAVE* HRAVE* H* HWALL* K* KINC* KMAX t LOC* L 

1 . M, NCHOKE • NERR* NFILE* NK* NMESH* NNN * NTAPE* PHI* PHIW 

1 » P* QR* 0* RCAP * RHOAV * RHO* RHOU* RROU* R* RUH 

1 * SIGMA* THETA* TRCL * U* VISC* W* WW * Z 

DIMENSION D I AM ( 2000 ) t P ( 2000 ) • E ( 2000 ) » HWALL(2000) 

1 * PHI ( 100 ) * S I GMA ( 100 ) * RCAP ( 100 ) * VI SC ( 100 ) • RHO ( 100 ) 

1 . R ( 100 ) * RHOU ( 100 ) * RROU ( 100 ) 

1 ♦ H ( 2*100 ) ♦ U ( 2 * 100 ) 

1 . RUPREV ( 100 ) ♦ RCPl(lOO)* RCP2 ( 100 ) • RCP3 ( 100 ) * RCP4(100) 

CALL CRISIS ( ZZ *RUPREV) 



CALCULATE THE VELOCTIY AT THE NEXT AXIAL STATION 
NMESHP = NMESH + 1 
IF(NNN) 10*10*20 
10 CONTINUE 
DRUT * 0.0 
U ( L * NMESHP) = 0.0 

AR = D I AM ( K-l ) *D I AM ( K“1 ) / ( DI AM ( K )#D I AM t K) ) 
DO 30 J=2.*NME5H 

RUPREV ( J ) = RROU(J) 

RROU(J) = RHOU(J) 


CORRECTION FOR RADIAL CONVECTION 
DRUTP = DRUT *R ( J-l ) /R ( J ) 

DRUT h RROU(J) - RUPREV ( J ) *AR + DRUTP 


62 

I F ( DRUT ) 

64*68 *66 

64 

GB 

= U ( L * J ) 


GO 

TO 68 

66 

GB 

= U ( L * J+l ) 

68 

I F( DRUTP 

) 70*74*72 

70 

GA 

= U t L * J-l ) 


GO 

TO 74 

72 

GA 

= U ( L * J ) 

74 

CONTINUE 



RADCON = 

DRUT*GB - DRUTP*GA 


1 + 1.0*U(L*J)*(RUPREV( J)*AR - RROU(J)) 


SMOOTHING THE RADIAL CONVECTION 

RCP4 ( J ) o RCP3 ( J ) 
RCP3 ( J ) = RCP2 ( J ) 
RCP2 ( J ) = RCPl(J) 
RCPKJ) = RADCON 
I F ( K - A) 42*42*30 
CONTINUE 


40 

42 



RCPltJ) = 0.0 
RCP2 ( J ) = 0.0 
RCP3 ( J ) = 0.0 
RCP4(J) = 0.0 
RUPREV(J) = RROU(J) 

30 CONTINUE 

U(L.NMESHP) = 0*0 - U( L#NMESH) 

UIM.NMESHP) = 0.0 
20 CONTINUE 

DO 50 J=2*NMESH 

RDR » R{ J)*DR 

DA * 6.2832*RDR 

VL = 0.0 - ( 1 . 0/RDR ) #( 

1 ( (R ( J+l ) + R(J) )/2*0)*( (VISC( J+l) + VISC(J))/2.0>* 

2 ( (U(L*J+1) - U(L*JM/DR) 

3 - ( ( R ( J ) + R(J-l) ) /2»0)#( (VISC(J) + VISCU-m/2.0)* 

4 I (U(L*J) - U(L* J-l ) )/DR) ) 

RADCON = 0.25 1 RCP1 ( J ) + RCP2 ( J ) + RCP3(J) +RCP4I J) ) 
* U(M*J) = U(L*J) - (DP + DZ*VL - RADCON )*( RUPREV ( J ) / 

2 (RROUf J)*RROU( J) ) ) 

50 CONTINUE 

U(M»1 I » U(M«2) 

RETURN 

END 


P 
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OA CHT0708 VAL 'WATSON SUBROUTINE ENERGY FOR HT0720 

SUBROUTINE ENERGY 

COMMON AMPS * DIAM* DIA* DP* DRDR * DR* DW ♦ DZ * EPS* E 
1 , EX* EXX, FZO * HAVE. HRAVE* H. HWALL* K* KINC. KMAX * LOC, L 

1 . M* NCHOKE* NERR* NFILE* NK* NMESH * NNN * NTAPE * PHI* PHIW 

1 * P. OR* 0* RCAP * RHOAV * RHO* RHOU* RROU* R* RUH 

1 > SIGMA* THETA* TRCL * U* VI SC* W* WW* Z 

DIMENSION D I AM ( 2000 ) * P ( 2000 ) • E(2000). HWALL(2000) 

1 * PHI < 100) * SIGMA(IOO)* RCAP ( 100 ) * VISC(IOO)* RHO(IOO) 

1 * R ( 100 ) * RHOU(IOO). RROU ( 100 ) 

1 , H ( 2 * 100 ) * U ( 2*100) 

1 * RCPl ( 100 ) * RCP2 ( 100 ) * RCP3I100). RCP4(100> 

1 . TE(IOO) 

CALL CRISIS ( ZZ *RCP1 ) 

C 

C CALCULATE THE ENTHALPY AT THE NEXT AXIAL STATION 
DRUT = 0.0 

TE ( 2 ) = U(M*2)*U(M*2) /2.0 

AR = D I AM ( K-l ) *DI AM ( K-l ) / ( DI AM ( K ) *D I AM ( K ) ) 

DO AO J=2*NMESH 
FJ = J 

R ( J ) = ( FJ— 1 • 5 ) *DR 
RDR = R ( J ) *DR 
DA = 6.2832*RDR 

CL = 0.0 - KPHKJ+1) - 2.0*PHI(J) + PHI ( J-l )) /DRDR 
1 + (PHKJ+l) - PHI (J-l) )/(2.0*RDR> ) 

RL = RCAP(J) 

OH = (E(K-1)*E(K-1)*SIGMA( J) ) 

TE ( J+l ) = U(M*J+l)*U(M*J+l)/2.0 
TEP = UtL.J )*U(L*J ) / 2 • 0 

C 

C CORRECTION FOR RADIAL CONVECTION 
DRUTP *> DRUT*R( J-l)/R( J) 

DRUT = RHOU(J) - RROUI J)*AR + DRUTP 


62 

IF(DRUT) 64*68*66 


64 

GB = H(L» J) 
GO TO 68 

+ TE ( J ) 

66 

GB = H { L * J+l ) 

+ TE ( J+l ) 

68 

I F ( DRUTP ) 70*74.72 


70 

GA = H(L*J-1) 
GO TO 74 

+ TE(J-l) 

72 

GA = H ( L * J ) 

+ TE ( J ) 

74 

CONTINUE 



RADCON = DRUT*GB - 

DRUTP*GA 


2 + 1«0*H(L*J)*(RROU(J) *AR - RHOU(J)) 

C 

C SMOOTHING THE RADIAL CONVECTION 



I 


RCPM J ) = RCP3(J> 

RCP3U) = RCP2 ( J ) 

RCP2 ( J ) = RCPl(J) 

RCPl(J) * RADCON 
IF(K-4) 20.20*22 
20 CONTINUE 

RCPl(J) = 0.0 
RCP2(J) = 0.0 
RCP3IJ) = 0.0 
RCP4U) s 0.0 
RROU(J) = RHOU(J) 

22 CONTINUE 

RADCON = 0 « 25 ( RCP1 ( J ) + RCP2(J) + RCP3IJ) +RCP4U)) 
H (M * J) = H ( L * J ) + (DZ*|OH-CL-RL) + RADCON )* (RROU( J )/ 

1 (RHOU(J)*RHOU(J) ) ) +TEP - TE ( J ) 

40 CONTINUE 

H(M.l) = H(M»2) 

RETURN 

END 


ON 

UJ 
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CHT0709 VAL WATSON SUBROUTINE OUTPT FOR HT0720 
^ SUBROUTINE OUTPT 

COMMON AMPS ♦ DIAM. DIA* DP* DRDR * DR* DW* DZ* EPS* E 
1 * EX* EXX* FZO * HAVE* HRAVE* H* HWALL* K» KINC* KMAX * LOC, L 

1 * M* NCHOKE * NERR* NFILE, NK* NMESH * NNN ♦ NTAPE * PHI* PHIW 

1 * P* OR* Q* RCAP • RHOAV * RHO* RHOU* RROU* R* RUH 

1 * SIGMA* THETA* TRCL * U* VISC* W * WW» Z 

DIMENSION D I AM ( 2000 ) * P < 2000 ) • E ( 2000 ) » HWAL.U 2000) 

1 » PHI ( 100 ) * S I GMA (100)* RCAP ( 100 ) * VISC(IOO)* RHO ( 100 ) 

1 * R ( 100 ) • RHOU ( 100 ) » RROU(IOO) 

1 * H ( 2 * 1 00 ) * U ( 2 * 100 ) 

CALL CRISIS ( ZZ.AA ) 

QT = Q+OR 
RUHA = RUH/W 
PRAD = QR/QT *100.0 
PRES = P(K) /1.013E5 

WRITE TAPE NTAPE* K« DIAM(K)* Z. AMPS* EDO* W»QT tHWALL DO 

1 .PRAD* NMESH* PRES * FZO *LOC *EPS *DW ♦ DP* DZ 

2 . HAVE* RUHA*HRAVE 

WRITE TAPE NTAPE* (R(J). H(M.J)* UIM»J>* RHOU(J)* J=1*NMESH) 

END FILE NTAPE 

WRITE OUTPUT TAPE 6,300*K* D I AM ( K ) *Z * AMPS * E ( K ) *W *QT . TRCL 

WRITE OUTPUT TAPE 6 * 299 * PRAD *HWALL ( K ) *PRE S* NMESH * LOC*FZO *DW* EPS * 

1 NFILE* EX 

WRITE OUTPUT TAPE 6* 298, HAVE, RUHA* HRAVE 

WRITE OUTPUT TAPE 6*301*(R(J), H(M*J!* U(M,J), RHOU(J), 

1 RIJ+25)* H ( M , J+25 ) * U(M*J+25)* RHOU1J+251* J=l*25 ) 

300 FORMAT ( 1H1 * 

1 15HAXIAL STATION = *14 * 16X* 

3 50H DC THERMAL ARC WITH AXIAL GAS FLOW , 

2 15HDIAMETER = •* E10.3* 10H METERS / 

1 IX* 15HAXIAL DIST = * E10.3* 10H METERS * 

2 50X* 15HCURRENT = , E10.3* 10H AMPS / 

1 IX* 15HVOLTAGE GRAD = ♦ E10.3* 10H VOLTS/M , 

2 50X* 15HFLOW RATE = , E10.3* 10H KG/SEC / 

2 IX, 15HWALL HEAT FLUX= * E10.3* 22H WATTS/M **2 * 

2 38X * 15HTRANS COOLING = * E10.3* 12H KG/SEC-M**2 ) 

299 FORMAT ( 1H * 

2 15HRADI ATION LOSS* * F10.4* 10H PERCENT * 

2 50X* 15HWALL ENTHALPY = * E10.3, 10H JOULES/KG / 

2 IX* 15HPRESSURE = , F10.4* 10H ATMOS * 

2 50X* 15HNMESH = , 14 / 

2 IX* 15HLOC = * 14* 16X * 

2 50X. 15HFZO = * E10.3 / 

2 IX, 15HDW = * E10 • 3 * 10X 

2 50X* 15HEPS = , E10.3 / 

2 IX* 15HFILE = * 1 4 *66X , 1 5HEX = *E10.3> 



298 FORMAT) 1H0* 

1 25HSPACE AVERAGE ENTHALPY = » E15.5* 10H JOULES/KG / 

2 IX* 25HMASS AVERAGE ENTHALPY * , E15.5* 10H JOULES/KG / 

3 IX* 25HAVERAGE ENERGY DENSITY = * E15.5* 12H JOULES/M**3 ) 

301 FORMAT ( 1H0* 3X*10H RADIUS *5X*10H ENTHALPY * 

1 5X ♦ 10H VELOCITY * 5X* 10H MASS FLUX * 

2 5X * 10H RADIUS *5X*10H ENTHALPY * 

3 5X * 10H VELOCITY * 5X* 10H MASS FLUX / 

4 IX, 3X ♦ 10H METERS *5X*10H JOULES/KG * 

5 5X * 10H M/S * 5X * 10H KG/S M**2 , 

4 5X*10H METERS *5X*10H JOULES/KG , 

5 5X * 10H M/S * 5X » 10H KG/S M**2 // 

9 ( 8 ( E15 *5 ) ) ) 

NFILE = NFILE + 1 

RETURN 

END 




I I Hill Illllll II I 
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APPENDIX C 


FORTRAN PROGRAM FOR THE SOLUTION OF THE ASYMMETRIC CONSTRICTED ARC 

WITH AN AXIAL FLOW OF GAS 


This appendix contains the Fortran II programs for the numerical solution 
of the asymmetric constricted arc with an axial flow of gas. The data 
required for input into these programs and a legend for the Fortran variables 
are also included. 

The function of each subroutine is described below. 


BO UNDO provides for the specification of the boundary conditions for the 

main program 

STATEP evaluates the state properties (except density) for all of the 
mesh points at each axial station 

WDOT evaluates the density for all mesh points at each axial station 

and calculates the mass flow rate by integration of the mass 
flux over the constrictor cross-sectional area 

ENERGY calculates the enthalpy for all mesh points at the next axial 
station from the energy equation 

MOM calculates the velocity for all mesh points at the next axial 

station from the momentum equation 

ITER provides the iteration of pressure to obtain the correct mass flow 

in the subsonic portion of the nozzle 

OUTPT prints and writes on magnetic tape the results of the main program 


A second main program entitled herein "Supersonic Continuation" uses the 
same subroutines to continue into the supersonic region of the nozzle with the 
exception that the following subroutine replaces ITER. 

ITERS provides the iteration of pressure to obtain the correct mass flow 

in the supersonic portion of the nozzle 
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.L IIIIIIIIIUI.il, 


INPUT DATA REQUIRED 


Card number 

Fortran variables 

Format 

1 

NMESH, NFILE, NT APE, NK 

4l4 

2 

CDIA 

E1C 

).3 

3 

AMPS 



4 

P(l) 



5 

EXX 



6 

CHWALL 



7 

FZO 



8 

EPS 



9 

EX 

\ 

/ 

10 

NCD 

14 

11 through 11 + NCD 

J, I, H(l, J, I) 

214, E10. 


These Fortran variables are described in the following legend. under the 
variables common to all programs or under the variables for program HT 0751* 
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INPUT DATA REQUIRED FOR SUPERSONIC CONTINUATION 


Card number 

Fortran variables 

Format 

1 

• WFILE(l ) , IWAPE(l) 

2l4 

2 

HFILE, NTAPE, KMAX, KINC 

414 

3 

EX 

E10.3 

k 

zc 

E10.3 

5 

THETA 

E10.3 


VARIABLES USED IN THE PROGRAMS FOR THE NUMERICAL SOLUTIONS OF THE 
ASYMMETRIC CONSTRICTED THERMAL ARC WITH AN AXIAL FLOW OF GAS 


Variable name 
AMPS 

DIAM 

DIA 

DP 

DW 

DZ 

EPS 

E 

EX 

EXX 

FZO 

H 

HWALL 

K 


Variables Common to all Programs 

Descr iption 

total current carried by the constrictor , A 

diameter of the constrictor (an array) , m 

diameter of the constrictor, m 

pressure drop between axial stations, N/m 2 

discrepancy between the mass flow rate at this axial 
station and the initial mass flow rate, kg/s 

incremental axial distance, m 

maximum allowable relative discrepancy of the mass flow 
rate 

axial voltage gradient, V/m 

factor by which the axial incremental distance is increased 
for the next axial station 

mass flow rate at the initial axial station, kg/s 

length of the first axial increment divided by the 
characteristic arc length, z G 

local enthalpy, j/kg 

enthalpy of the gas at the constrictor wall, J/kg 
axial station number 
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Variable name 


Description 


LOC 

L 

M 

NCHOKE 

NFILE 

NK 

HMESH 

HM 

HTAPE 

PHI 

PHIW 

P 

QR 

Q 

RCAP 

RHOAV 

RHO 

RHOU 

SIGMA. 

U 


number of iterations required to satisfy the continuity 
equation (i.e., the number of iterations required such 
that the flow rate at this axial station is sufficiently 
near the initial flow rate) 

relative axial station number 

second relative axial station number 

flag indicating that the flow is choked 

file number on the magnetic tape used to store the solutions 

maximum number of axial stations to be calculated 

number of radial increments from the center line to the 
constrictor wall 

flag indicating that this is the first iteration on pres- 
sure drop (zero indicates first iteration, 1 indicates 
all iterations thereafter) 

magnetic tape on which the solution is stored 

local thermal conductivity potential, W/m 

thermal conductivity potential of the gas at the constrictor 
wall, W/m 

static pressure at this axial station, N/m 2 

heat transfer rate to the constrictor wall that is due to 
radiation, w/m 2 

heat transfer rate to the constrictor wall that is due to 
thermal conduction, W/m 2 

local radiation, W/m 3 

space average density, kg/m 3 

local density, kg/m 3 

local product of density and velocity at this station, 
kg/sm 2 

local electrical conductivity, l/R-m 
local velocity, m/s 


TO 


Variable name 


Descri p tion 


vise 

W 

ww 

z 

CWF 

DZMAX 

RHOA 

ZO 

CDIA 

CHWALL 

FJ 

FNI 

I 

J 

NCD 

NI 


local viscosity, Ns/m 2 

mass flow rate at the initial axial station, kg/s 
mass flow rate at this axial station, kg/s 

axial distance of this axial station from the initial axial 
station, m 

Variables Local to Program HT-0750 

scaling factor that changes the magnitude of the velocity 
profile to obtain the desired initial flow rate 

upper limit of the axial incremental distance (the incre- 
mental axial distance is never allowed to exceed this 
value in order to keep the solution stable) , m 

space average of the density at this axial station, kg/m 3 

characteristic arc length, z Q 

Variables Local to Program HT-O 75 I 
constrictor wall diameter, m 

enthalpy of the gas at the constrictor wall, j/kg 
floating index for the radial position 

number of increments in the azimuthal direction (floating) 
index for the azimuthal position 
index for the radial position 

number of data cards used to specify the initial enthalpy 
distribution 

number of azimuthal increments 


Variables Local to Program HT-0752 
DA incremental cross-sectional area, m 2 

DRDR square of the incremental radial distance, DR, m 2 
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Variable name 

m 

USE 

FJ 

FNI 

I 

J 

NERR 

NI 

QRR 

R 

RW 

SS 

SW 


Des cription 

incremental radial distance, m 

incremental azimuthal distance, radians 

floating index for the radial position 

number of increments in the azimuthal direction (floating) 

index for the azimuthal position 

index for the radial position 

error flag for the gas property subroutines 

number of azimuthal increments 

running total of the product of radiation and area, W/m 
radial position, m 
dummy variable not used 

running total of the product of the electrical conductivity 
and area, m/Q. 

dummy variable not used 


DA 

DRDR 

m 

DTH 

FJ 

FNI 

I 

J 

NERR 

RI 

R 


Variables Local to Programs HI -0753 and 0754 
incremental cross-sectional area, m 2 
square of the incremental radial distance, DR, m 2 
incremental radial distance, m 
incremental azimuthal distance, radians 
floating index for the radial position 

number of increments in the azimuthal direction (floating) 

index for the azimuthal position 

index for the radial position 

error flag for the gas property subroutines 

number of azimuthal increments 

radial position, m 
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Variable name 

CL 

DA 

DR DR 

DR 

DTH 

DZRU 

FJ 

FNI 

ICC 

IC 

IMM 

IM 

IPP 

IP 

I 

J 

NI 

OH 

PHIMM 

PHIM 

PHIO 


Variables Local to Program HT-0755 

Description 

conduction losses from this incremental volume, W/m 3 

incremental cross-sectional area, m 2 

square of the incremental radial distance, DR, m 2 

incremental radial distance, m 

Incremental azimuthal distance, radians 

local DZ divided by the local RHOU, m 3 s/kg 

index for the radial incremental distance (floating) 

number of azimuthal increments 

index for changing the number of azimuthal increments 

second index for changing the number of azimuthal 
increments 

index for the previous azimuthal position 

index for the azimuthal position at the previous radial 
position 

index for the next azimuthal position 

index for the azimuthal position at the next radial 
position 

index for the azimuthal position 

index for the radial position 

number of azimuthal increments 

ohmic heating within this volume element, W/m 3 

thermal conductivity potential at the previous radial 
position, W/m 

thermal conductivity potential at the previous azimuthal 
position, W/m 

thermal conductivity potential at this radial and 
azimuthal position, W/m 
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Variable nam e 
PHIPP 

PHIP 

RDR 

RDTH 2 

RL 

R 

SLOPE 

DA 

DRDR 

DR 

DTK 

DZRU 

PJ 

FNI 

ICC 

IC 

IMM 

IM 

IPP 

IP 


Descri ptio n _ _ 

thermal conductivity potential at the next azimuthal 
position, W/m 

thermal conductivity potential at the next radial 
position, W/m 

product of the local radius and incremental radial distance, 
m 2 

square of the product of radius and incremental azimuthal 
distance, m 2 

radiation losses from this incremental volume, w/m 3 
radial position, m 

gradient of the thermal conductivity potential at the wall, 
W/m 2 

Variables Local to Program HT-O756 
incremental cross-sectional area, m 2 
square of the incremental radial distance, DR, m 2 
incremental radial distance, m 
incremental azimuthal distance, radians 
local DZ divided by the local RHOU, m 3 s/kg 
index for the radial incremental distance (floating) 
number of azimuthal increments 

index for changing the number of azimuthal increments 

second index for changing the number of azimuthal 
increments 

index for the previous azimuthal position 

index for the azimuthal position at the previous radial 
position 

index for the next azimuthal position 

index for the azimuthal position at the next radial 
position 
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Variable name 


I 

J 

NI 

RDR 

RDTH2 

R 

UMM 

UM 

UO 

UPP 

UP 

VL 

VMM 

VM 

VO 

VPP 

VP 


DDP 

IND 

N 


Description 

index for the azimuthal position 

index for the radial position 

number of azimuthal increments 

product of the local radius and incremental radial 
distance , m 2 

square of the product of radius and incremental azimuthal 
distance , m 2 

radial position, m 

velocity at the previous azimuthal position, m/s* 
velocity at the previous radial position, m/s 
velocity at this radial and azimuthal position, m/s 
velocity at the next azimuthal position, m/s 
velocity at the next radial position, m/s 
viscous losses from this volume increment, N/m 3 
viscosity at the previous azimuthal position. Ns/m 2 
viscosity at the previous radial position. Ns/m 2 
viscosity at this radial and azimuthal position. Ns/m 2 
viscosity at the next azimuthal position. Ns/m 2 
viscosity at the next radial position, Ns/m 2 

Variables Local to Program HT-O757 

amount that the pressure drop has changed for each 
iteration at this axial station, N/m 2 

flag used in adjusting the pressure drop during the 
iterations at this station 

index for the iterations 
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Variables Local to Program HT-O758 


Variable 

FJ 

FNI 

J 

NI 

PRAD 

qr 

KINC 

KMAX 

NFILE 1 

NI 

NN 

NT APE 1 
PX 

qx 

THEIAR 

THETA 

zc 

DDP 


name Description 

index for the radial position (floating) 
number of azimuthal increments (floating) 
index for the radial position 
number of azimuthal increments 

percentage of the constrictor wall heat flux that is due to 
radiation 

total heat transfer rate to the constrictor wall, W/m 2 

Variables Local to Program HT-O78O 

number of axial stations between printout of results 

maximum number of axial stations to be calculated 

magnetic tape file number indicating the last file of the 
subsonic results 

number of azimuthal increments 

index for printing results 

magnetic tape on which the subsonic results are stored 

pressure at the last axial station of the subsonic flow 
regime, N/m 2 

heat flux to the constrictor wall at the last axial station 
of the subsonic regime, W/m 2 

half-angle of divergence of the supersonic nozzle, radians 
half-angle of divergence of the supersonic nozzle, deg 
axial position at choking, m 

Variables Local to Program HT-O781 

amount that pressure drop has changed for each iteration 
at this axial station, N/m 2 
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Variable 

USD 

N 


name Description 

flag used in adjusting the pressure drop during the 
iterations at this station 

index for the iterations 
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-J CHT0750 VAL WATSON DC ARC - AXIAL FLOW - REAL GAS 

^ COMMON AMPS » W, WW* DP* DZ* EPS* FZO* Z* Q* QR * DW* DIA* EX* EXX* 

1 H* PHI* SIGMA* RCAP ♦ RHO* VISC* HWALL* RHOAV * RHOU* PHIW* 

1 P* E* DIAM* U* 

1 NMESH * K » L* M* LOC* NFILE* NTAPE ♦ NCHOKE • NK* NNN 

DIMENSION DIAM(IOOO)* P(IOOO)* E(IOOO)* HWALL(IOOO) 

DIMENSION H ( 2 » 1 2 *48 ) * U(2*12*481* PHI(12*48)* S IGMA ( 12 *48 ) * 

1 RCAP (12*48). V I SC (12 *48 ) * RHO (12*48)* RHOU(12*48) 

CALL CRISIS ( ZZ * AAA ♦ AMPS * NNN) 

CALL BOUNDC 
K = 1 

DIA = DIAM(K) 

L = 1 
M = 1 
Z = 0.0 
CALL STATEP 
CALL WDOT 

CWF = EXX/WW 
DO 20 J=1 *NMESH 

DO 18 1=1*48 

18 U ( 1 ♦ J * I ) = CWF*U( 1*J* I ) 

20 CONTINUE 

CALL WDOT 
W = WW 

ZO = (W* 3111.0) /3.1416 
DZ = FZO#ZO 
CALL RHOAVE 
RHOA = RHOAV 

CALL LOCATE (NFILE.NTAPEI 
CALL OUTPT 
M = 2 
K = 2 

DIA = DIAM(K) 

DP = 0.0 

P(K) = P(K-l) + DP 
CALL ENERGY 
CALL RHOAVE 

DP =W*W* ( RHOAV- RHO A ) / ( ( ( RHOA#3 • 1416*D I A*DI A ) /4.0 ) **2 ) 

CALL ITER 
Z = Z + DZ 
CALL STATEP 
CALL OUTPT 
DO 6 K = 3 * NK 
DIA = DIAM(K) 

DZMAX =• 0.2*DIA*DIA#RHOU(3*1 >*200.0/ ( FLOATFf NMESH )*FLOATF (NMESH) ) 
IF(DZ-DZMAX) 40*42.42 
40 DZ = DZ*FX 


I 


42 CONTINUE 

M = 1 
L “ L + 1 
IF(3-L) 1*1*2 

1 M = 2 
L = 1 

2 CONTINUE 
CALL ENERGY 
CALL ITER 

I F ( NCHOKE ) 4*4*3 

3 WRITE OUTPUT TAPE 6* 202# K* DW#U(M*1*1) 

GO TO 8 

4 CONTINUE 
CALL STATEP 
Z = Z + DZ 
CALL OUTPT 

6 CONTINUE 
8 CONTINUE 
REWIND NTAPE 
REWIND 8 

202 FORMAT ( 1 HO * 18HFLOW CHOKED AT K = # 14* 10X# 

1 20HFLOW RATE ERROR IS * 2PF12«7# 10H PERCENT 

2 14 HCL VELOCITY = *0PF10#1* 8H M/SEC ) 

CALL EXIT 

END 



VO 




CO CHT0751 VAL WATSON SUBROUTINE BOUNDC FOR HT0750 

° SUBROUTINE BOUNDC 

COMMON AMPS » W> WW* DP* DZ* EPS* FZO* Z* Q* OR* DW* DIA* EX* EXX* 

1 H* PHI* SIGMA* RCAP ♦ RHO* VISC* HWALL* RHOAV* RHOU* PHIW* 

1 P* E* DIAM* U. 

1 NMESH * K » L* M* LOC* NFILE. NT APE ♦ NCHOKE * NIC* NNN 

DIMENSION DIAM! 1000) ♦ P ( 1000 ) * E<1000)* HWALL(IOOO) 

DIMENSION H ( 2 * 12 *48 ) * U(2*12»48)» PHI(12*48)* SIGMA ( 12 *48 ) * 

1 RCAP ( 1 2 * 48 ) • V ISC ( 1 2 *48 ) » RHO(12*48l* RH0U(12*48) 

CALL CRISIS ( ZZ *AAA) 

READ INPUT TAPE 5* 99* NMESH* NFILE* NTAPE * NK 
99 FORMAT (414) 

READ INPUT TAPE 5* 100* CDIA* AMPS* P(l)* EXX* CHWALL* FZO* EPS*EX 
100 FORMAT ( E10* 3 ) 

Pll) = P ( 1 ) *1 .013EK 
DO 5 K = l.NK 

DIAM(K) a CDIA 
HWALL ( K ) = CHWALL 

5 CONTINUE 
N I = 6 

DO 10 J = 1 *NMESH 
DO 6 I = 1 » N I 

H(1*J*I) = HWALL ( 1 ) 

6 U ( 1 * J . I ) = 100.0 
FJ = J 

FNI = N I 

IF ( (FNI/FJ1-4.0) 8*8*10 
8 N I = N 1*2 

10 CONTINUE 

READ INPUT TAPE 5* 102* NCD * (J* I. H(1*J*I)* K=1*NCD ) 

102 FORMAT! I4/< 214. E10. 3) ) 

RETURN 

END 
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CHT0752 VAL .WATSON SUBROUTINE STATEP FOR HT0750 
SUBROUTINE STATEP 

COMMON AMPS* W* WW* DP* DZ* EPS* FZO* Z* O* OR* DW* DIA* EX* EXX* 
1 H* PHI* SIGMA* RCAP * RHO* VISC* HWALL* RHOAV ♦ RHOU* PHIW* 

1 P, E* DIAM* U* 

1 NMESH. K* L* M* LOC* NFILE* NT APE * NCHOKE * NK* NNN 

DIMENSION DIAMt 1000) * P(IOOO)* E(IOOO)* HWALL(IOOO) 

DIMENSION H ( 2 ♦ 12 *48 ) * U(2*12*48)* PHI(12*48)* S IGMA ( 1 2 *48 ) * 

1 RCAP (12*48) ♦ VI SC ( 1 2 *48 ) • RHO ( 1 2 *48 ) • RH0U(12*48) 

CALL CR I SI S ( ZZ * AAA ) 

OR = QRR 

PRES=P(K)/1.013E5 

CALL NTAB ( PRES ♦ HWALL ( K ) * PHIW* SW* RW* VISC ( NMESH+1 *11* 8* NERR) 
DR = DIA/(2.0*FLOATF(NMESH 1-1.0) 

DRDR = OR*DR 
R a o.O 
J = 1 
I = 1 

CALL NT AB ( PRES* H(M*J*D* PHI(J*D* SIGMA ( J* I ) * 

1 RCAP ( J * I ) * VISC(J*I1* 8* NERR) 

PHI ( 1 * 2 ) = PHI ( 1 *1 ) 

PHI ( 1*3) a PHin.l) 

VISC (1*2) = VISC ( 1 *1 ) 

V ISC ( 1 *3 ) = VISCt 1*1) 

SS = SIGMAt J* I ) *3. 1416*DRDR#0.25 
ORRs RCAP(J*I )*3.1416*DRDR*0.25 
NI » 6 

DO 14 J a 2*NMESH 
R = R + DR 

DTH a 6.2832/FLOATF (NI ) 

DA = R*DR*DTH 

DO 10 I = 1 *NI 

CALL NTAB(P( K) * H(M*J*I)> PHI(J*D* SIGMA(J*D* 

1 RCAP(J*D* VISC t J * I ) * 8* NERR) 

SS » SS + SIGMA ( J ♦ I ) *DA 
QRRa QRR+ RCAP ( J * I ) *DA 
10 CONTINUE 

FJ = J 
FNI = NI 

IF M FNI /FJ ) - 4.0) 12*12*14 

12 NI a N 1*2 

14 CONTINUE 

E(K) = AMPS/SS 
QRR=QRR/(3.1416*DIA) 

I F ( NERR ) 22*22*20 
20 WRITE OUTPUT TAPE 6* 110* K 
22 CONTINUE 



110 FORMAT ( 1H0 * 24HEXCEEDED TABLES AT K = ♦ 14) 

RETURN 
END 


CHT0753 VAL WATSON SUBROUTINE WDOT FOR HTO750 
SUBROUTINE WDOT 

COMMON AMPS* W* WW* DP* DZ* EPS* FZO* Z* Q* QR * DW* DIA* EX* EXX* 
1 H, PHI. SIGMA. RCAP * RHO* VISC* HWALL* RHOAV. RHOU, PHIW* 

1 P*E*DIAM*U* 

1 NMESH* K* L* M* LOC* NFILE* NTAPE * NCHOKE* NIC. NNN 

DIMENSION D I AM ( 1000 ) * P(IOOO)* E(lOOO)* HWALL(lOOO) 

DIMENSION H(2*12*48). U ( 2 . 12*48). PHI { 12*48 ) • S IGMA ( 1 2 *48 ) . 

1 RCAP (12*48) » VISC ( 12 *48) * RH0(12*48)* RH0U(12*48) 

CALL CRI 5 IS ( ZZ ♦ AAA ) 

DR = DIA/(2»0#FLOATF(NMESH)-1.0) 

DRDR = DR*DR 
R = 0.0 
J = 1 
I = 1 

PRES=P(K)/1.013E5 

CALL NRHO ( PRES * H(M*J.I). RHO ( J* I ) * 8* NERR) 

RHOU ( J * I ) = RHO ( J * I ) *U ( M* J * I ) 

WW □ RHOU ( J * I ) *3 . 1416*DRDR*0« 25 
N I = 6 

DO 14 J = 2 *NMESH 

p — p ^ pp 

DTH = 6. 2832/FLO ATF (NI ) 

DA = R*DR*DTH 
DO 10 I * 1 *NI 

CALL NRHO ( PRES * H(M*J.I)» RHO(J»I). 8. NERR) 

RHO.U ( J * I ) = RHO(J.I)*U(M»J.I) 

WW a WW + RHOU ( J * I ) *DA 
10 CONTINUE 

FJ = J 
E NI = NI 

IF ((FNI/FJ) - 4.0) 12*12*14 

12 NI = N 1*2 

14 CONTINUE 

I F ( NERR ) 22*22*20 
20 WRITE OUTPUT TAPE 6* 110* K 
22 CONTINUE 

110 FORMAT ( 1H0 * 24HEXCEEDED TABLES AT K = . 14) 

RETURN 

END 



CHT0754 VAL WATSON SUBOUTINE RHOAVE FOR HT0750 
SUBROUTINE RHOAVE 

COMMON AMPS , W» WW* DP» DZ* EPS, FZO* Z* Q» QR * DW, DIA* EX* EXX* 
1 H* PHI, SIGMA* RCAP* RHO* VISC* HWALL* RHOAV * RHOU* PHIW* 

1 P* E* DIAM* U. 

1 NMESH , K » L* M, LOC* NFILE, NTAPE * NCHOKE • NK* NNN 

DIMENSION D I AMI 1000 ) * P(IOOO), E ( 1 000 ) » HWALLI1000) 

DIMENSION H ( 2 * 12 ,48 ) • U ( 2 * 12*481* PHI(12*481* S IGMA ( 12 *48 ) * 

1 RCAP(12,48), VI SC (12*48) * RH0(12*48)i RH0U(12*48) 

CALL CR I S IS ( ZZ » AAA ) 

DR = DIA/ (2 «0*FL0ATF (NMESHl-1,0) 

DRDR = DR*DR 
R = 0.0 
J = 1 
I = 1 

PRES=P (KJ/1.013E5 

CALL NRHO ( PRES * H(M*JtI), RHO(J*I)* 8, NERR) 

RHOAV = 3.1416*DRDR*0»25*RHO(J»I ) 

N1 = 6 

DO 14 J = 2 *NMESH 
R = R + DR 

DTH = 6.2832/FLOATF ( N I ) 

DA = R*DR*DTH 
DO 10 I = 1 , N I 

CALL NRHO ( PRES * H(M,J*I), RHOIJ, I)* 8, NERR) 

RHOAV = RHOAV +DA*RHO{J*I) 

10 CONTINUE 

FJ = J 
FNI = N I 

IF ((FNI/FJ) - 4.0) 12*12*14 

12 N I = N 1*2 

14 CONTINUE 

RHOAV = RH0AV/(3,1416*DIA*DIA*0.25) 

I Ft NERR ) 22 *22*20 
20 WRITE OUTPUT TAPE 6* 110* K 
22 CONTINUE 

110 FORMAT ( 1 HO » 24HEXCEEDED TABLES AT K = , 14) 

RETURN 

END 




1 

a> CHT0755 VAL WATSON SUBROUTINE ENERGY FOR HT0750 

SUBROUTINE ENERGY 

COMMON AMPS , W» WW* DP* DZ» EPS* FZO * Z* Q* OR* DW* DIA* EX* EXX* 

1 H* PHI, SIGMA, RCAP * RHO, VISC* HWALL* RHOAV * RHOU* PHIW* 

1 P* E* DIAM* U* 

1 NMESH , K* L, M* LOC* NFILE, NTAPE » NCHOK.E ♦ NK* NNN 

DIMENSION DIAMt 1000) * P(1000)* E(IOOO)* HWALL (1000) 

DIMENSION H ( 2 , 12 *48 ) , U(2, 12,48), PHI ( 12,4-8) • SIGMA ( 12 ,48 ) * 

1 RCAP (12*48) , V ISC ( 12 *48 ) * RH0(12*48)* RH0U(12*48) 

CALL CRI SIS ( ZZ , AAA ) 

0 = 0,0 
I * 1 
J = 1 

DR * D I A / ( 2 »0*FLOATF ( NMESH) -1.0 ) 

DRDR * DR*DR 
R = 0.0 

DZRU = DZ/RHOU ( J * I ) 

CL n ( 6»0*PHI (1*1) - PHI { 2 * 1 ) - PHI(2*2) - PHI(2*3) - PHI (2*4) 

1 - PHI ( 2 ♦ 5 ) - PHI ( 2 *6 ) ) / ( 6.0*0.25 *DRDR ) 

RL = RCAP ( 1*1) 

OH = (E(K-1)*E(K-1)*SIGMA(J,I) ) 

HIM, J*I ) = H ( L* J * I > + (DP/RHO{ J* I ) ) 

1 +DZRU* ( OH -CL -RL) 

ICC = 1 
NI □ 6 

DO 40 J = 2, NMESH 
FJ □ J 
FNI = NI 

IF ((FNI /FJ )— 4. 0 ) 1*1*3 

1 IC = 1 

3 CONTINUE 

R = R + DR 
RDR = R*DR 

DTH = 6.2832/FLOATF(NI ) 

DA = R*DR*DTH 
RDTH2 = R*R*DTH*DTH 
DO 30 I = 1 ,NI 

IP = I 

IM = I 

IPP =1+1 
I MM = I - 1 
IF(I-I) 2*2*4 

2 IMM = NI 

4 I F( NI - I) 6*6*8 

6 IPP = 1 

8 IF(IC) 12,12,10 

10 IP = 2*1 



12 I F ( ICC ) 16*16*14 

14 IM = (I+1I/2 

16 CONTINUE 

DZRU = DZ/RHOU ( J * I ) 

PHIO = PH I ( J * I ) 

PHIP = PH I ( J + l » I P ) 

PHIM a PHI ( J-l * IM ) 

PHIPP a PHI ( J * I PP ) 

PHIMM a PHI ( J * IMM ) 

I F ( NMESH-J ) 20*20.22 

20 PHIP = 2.0*PHlW - PH I ( J * I ) 

SLOPE a ( PH I W -PHIP ) *2 »0/DR 

IF (Q-SLOPE) 21,22*22 

21 Q = SLOPE 

22 CONTINUE 

CL = 0.0 - ( (PHIP-PHIM)/(2.0*RDR) 

1 + ( PHIP-2»0*PH IO+PHIM) /DRDR 

1 + (PHIPP-2.0*PHI0+PHIMM)/RDTH2) 

RL b RCAP ( J * I ) 

OH = ( E ( K-l ) *E ( K-l ) *S IGMA ( J * I ) ) 

H (M, J , I ) = H ( L , J » I ) + (DP/RHOt J*I ) ) 

1 + DZRU* ( OH - CL - RL) 

30 CONTINUE 

ICC = 0 

IF(IC) 40*40*32 
32 ICC a 1 

IC a o 
N I = 2*NI 

40 CONTINUE 
RETURN 
END 
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00 CHT0756 VAL WATSON SUBROUTINE MOM FOR HT0750 

SUBROUTINE MOM 

COMMON AMPS. W. WW. DP. DZ » EPS. FZO. Z. Q, QR * DWt DIA. EX. EXX. 

1 H, PHI, SIGMA. RCAP » RHO. VISC. HWALL. RHOAV. RHOU , PHIW. 

1 P, E » DIAM, U» 

1 NMESH, K* L. M* LOC, NFILE. NT APE ♦ NCHOKE . NK* NNN 

DIMENSION D I AM ( 1000 ) * P(IOOO)* E(IOOO), HWALL(IOOO) 

DIMENSION H ( 2 , 12 ,48 ) » U<2, 12.48). PHK12.48). S IGMA ( 12 ,48 ) , 

1 RCAP (12.48) » VISCI12 .48) , RH0(12,48), RHOU( 12 *48 ) ,RROU( 12 *48 > 

CALL CRISIS (ZZ.AAA) 

IF (NNN) 99,99*102 
99 IF (NMESH - 7) 200.201,201 

200 N I = 24 
GO TO 202 

201 N I = 48 

202 CONTINUE 

DO 101 J = 1 * NMESH 

DO 100 1 = 1 * N I 

RROU ( J * I ) = RHOU( J * I ) 

100 CONTINUE. 

101 CONTINUE 

102 CONTINUE 
I = 1 

J = 1 

DR = DIA/(2.0*FLOATF(NMESH>-1.0) 

DRDR = DR*DR 
R = 0.0 

DZRU = DZ/ RROU(J.I) 

VL = VISC(J»I)*I6.0*U(L*1»1) - U ( L * 2 • 1 ) - U(L*2*2) - U( L* 2 *3 ) - 
1 U(L,2»4) - UfL.2.5) -U ( L * 2 *6 )) / ( 6.0*3. 1416*DRDR ) 

U ( M , J » I ) = U(L.J.I) - ( DP/RROU ( J* I ) ) - DZRU*VL 

U ( M , 1 » 2 ) = U(M.l.l) 

U ( M » 1 ,3 ) = U ( M * 1 * 1 ) 

ICC = 1 
N I = 6 

DO 40 J = 2. NMESH 
FJ = J 
FNI a NI 

IFUFNI/FJ) -4.0) 1*1.3 

1 IC = 1 

3 CONTINUE 

R = R + DR 
RDR = R*DR 

DTH = 6.2832/FLOATF ( NI ) 

DA = R*DR*DTH 
RDTH2 = R*R*DTH*DTH 
DO 30 I = 1 »NI 


IP = I 

IM = I 

IPP = I + 1 

I MM = I - 1 

IF(I-l) 2*2*4 
2 I MM = NI 

4 I F ( N I - I) 6.6*8 

6 IPP b 1 

8 IF(IC) 12.12.10 

10 IP s 2*1 

12 I F ( ICC ) 16.16.14 

14 IM » (1+11/2 

16 CONTINUE 

DZRU = DZ/RROU ( J . I ) 

UO s U(L.J.I) 

UP = U ( L. J+l . I P 1 
UM s U(L.J-l.IM) 

UPP = U(L.J.IPP) 

UMM a U(L.J.IMM) 

VO b VISC(J.I) 

VP = VISC(J+1.IP) 

VM a VISC(J-l.IM) 

VPP = VISC(J.IPP) 

VMM = VISC(J.IMM) 

IF(NMESH-J) 20.20.22 

20 UP = 0.0 - U(L.J.I) 

VP a VISCtNMESH+l. 11*2.0 - VISC(J.I) 

22 CONTINUE 

VL = ( ( (R+DR/2.0)*( (VP+VO)/2.0)*(UO-UP) + 

1 ( R-DR/2 .0 ) *( ( VO+VM 1/2.01* (UO“UM ) ) / ( RDR*DR 1 

2 + ( ( (VPP+VO)/2.0)*(UO-UPP) + ( (VO+VMM)/2.0)*(UO-UMM) 1 

3 / ( RDTH2 1 1 

U(M.J.I) = U(L.J.I) - ( DP/RROU ( J. I ) 1 - DZRU*VL 

30 CONTINUE 

ICC = 0 

IF(IC) 40.40.32 
32 ICC * i 

IC = 0 
NI 8 2 *N I 

40 CONTINUE 
RETURN 
END 


A 


CD 

VO 



CHT0757 VAL WATSON SUBROUTINE ITER FOR HT0750 
SUBROUTINE ITER 

COMMON AMPS » W» WW* DP* DZ* EPS* FZO* Z* Q* QR* DW* DIA* EX* EXX* 
1 H. PHI. SIGMA. RCAP ♦ RHO* VISC* HWALL* RHOAV. RHOU* PHIW* 

1 P* E* DIAM* U* 

1 NMESH * K * L* M* LOC* NFILE* NTAPE * NCHOKE * NK. NNN 

DIMENSION D I AM ( 1 000 ) * P 1 1000 i • E ( 1000 ) • HWALL(IOOO) 

DIMENSION H ( 2 * 12 *48 ) * U(2*12*48)* PHK12.48)* S IGMA ( 12 *48 ) * 

1 RCAP (12*48) * VISC( 12 *48 ) » RHO (12*48)* RH0U(12*48) 

CALL CRI SIS t ZZ * AAA ) 

IND = 0 
NCHOKE = 0 
DDP » DP 

P ( K ) = P(K-l) + DP 
LOC = 0 
NNN » 0 

DO 4 N » 1*30 

LOC » N 
CALL MOM 
NNN a 1 

CALL WDOT 

DW a (WW - W)/W 

IF(ABSFIDW) - EPS) 10*20*20 


20 

IF 1 IND) 22 *22*1 


1 

DDP a DDP/2.0 


22 

IF(DW) 2*10*3 


2 

DP = DP + DDP 



P ( K ) a P(K-l) 
GO TO 4 

+ DP 

3 

DP a DP - DDP 



P ( K ) = PtK-1) 
IND a 1 

+ DP 

4 

CONTINUE 
NCHOKE a l 


10 

CONTINUE 

RETURN 

END 




I 
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^ CHT0758 VAL WATSON SUBROUTINE OUTPT FOR HT0750 

SUBROUTINE OUTPT 

COMMON AMPS. W. WW. DP. DZ. EPS. FZO* Z. Q. OR* DW* DIA* EX, EXX* 

1 H, PHI, SIGMA. RCAP » RHO, VISC* HWALL* RHOAV. RHOU* PHIW* 

1 P, E, DIAM, U, 

1 NMESH » K, L. M» LOC. NFILE* NT APE » NCHOKE * NK, NNN 

DIMENSION DIAM(IOOC), P(IOOO), E(IOOO)* HWALL(IOOQ) 

DIMENSION H ( 2 » 12 .43 ) . U<2»12»48), PH1(12*48), SIGMA (12 *48 ) « 

1 RCAP (12.48), VISC ( 12 .48 ) » RH0(12*48), RH0U(12*48) 

CALL CRISIS (ZZ.AAA) 

OT = Q+QR 

PRAD = QR/QT *100.0 

WRITE TAPE NTAPE, K, DI AM ( K. ) * Z, AMPS* E (K) * W »QT.*HWALL (K) 

1 .PRAD* NMESH, P I K) *FZO *LOC*E PS *DW* DP, DZ 

IF (NMESH - 7) 11,12,12 

11 N I = 24 
GO TO 13 

12 N I = 48 

13 CONTINUE 

WRITE TAPE NTAPE, NI, t(H(M.J,I)» U(M*J*I)* RHOU( J*I), I=1*NI), 

1 J = l, NMESH) 

END FILE NTAPE 

WRITE OUTPUT TAPE 6,300, K» DI AM ( K ) *Z ,AMPS*E ( K ) ,W *OT , HWALL ( K ) 

P(K) = P ( K) /1.013E5 

WRITE OUTPUT TAPE 6 , 299, PRAD *NMESH*P ( K )* FZOiLOC* EPS * DW*£X « NFI LE 

P ( K ) = P ( K ) *1 .013E5 

WRITE OUTPUT TAPE 6 » 30 1 *H (M , 1 , 1 ) 

NI = 6 

DO 2 J = 2, NMESH 

WRITE OUTPUT TAPE 6 *302 * (HI M * J » I ) *1 =1»NI) 

FJ = J 
FNI = NI 

IF ( (FNI/FJ) -4.0) 1,1,2 

1 NI = 2#N I 

2 CONTINUE 

WRITE OUTPUT TAPE 6*301* RHOU(l*l) 

NI = 6 

DO 4 J = 2, NMESH 

WRITE OUTPUT TAPE 6,302, fRHOU(J.I)* 1=1 *NI ) 

FJ = J 
FNI = NI 

I F ( (FNI/FJ) -4,0) 3*3*4 

3 NI = 2*NI 

4 CONTINUE 
300 FORMAT ( 1H1 , 

1 15HK = ,14 , 16X * 

3 50H DC THERMAL ARC WITH AXIAL GAS FLOW 



J 



2 15HDI AMETER 

1 IX* 15HAXIAL DIST 

2 50X » 15HCURRENT 

1 IX* 15HV0LTAGE GRAD = 

2 50X * 15HFL0W RATE 

2 IX* 15HWALL HEAT FLUX = 

2 38X » 15HWALL ENTHALPY = 

299 FORMAT ( 1H * 

2 15HRAD I AT ION LO SS = 

2 50X * 15HNMESH 

2 IX* 15HPRESSURE 

2 SOX * 15HFZO 

2 IX* 15HLOC 

2 SOX* 15HEPS 

2 IX* 15HDW 

2 50X * 15HEX 

2 IX* 15HFILE 

301 FORMAT ( 1H0* E12 • 3 ) 

302 FORMAT ( 1H0* (6E12.3/) ) 

NFILE = NFILE +1 
RETURN 
END 


E10.3* 10H METERS 

E10.3* 10H METERS 

E10.3* 10H AMPS 

E10 *3 * 10H VOLTS/M 

E10 1 3 * 10H KG/SEC 

E10.3* 22H W/M**2 
E10.3* 10H J/ KG 

F10.3* 10H PERCENT 

14 / 

E10* 3 * 10H ATMOS 

E10.3 / 

14. 16X * 

E10.3 / 

E10.3.10X 
E10.3 / 

14) 


i 

* 

/ 

* 

/ 

PREV STA 
) 




* 


VO 

U) 


CHT0780 VAL WATSON DC ARC - SUPERSONIC CONTINUATION - REAL GAS 

COMMON AMPS* W» WW. DP* DZ* EPS* FZO* Z* Q* QR* DW» DIA* EX* EXX. 
1 H* PHI* SIGMA* RCAP * RHO* VISC* HWALL* RHOAV * RHOU* PHIW* 

1 P. E* DIAM* U* 

1 NMESH * K * L* M* LOC* NFILE* NTAPE * NCHOKE * NIC* NNN 

DIMENSION D I AM ( 1000 ) * P(IOOO)* E(1000)* HWALL(IOOO) 

DIMENSION H ( 2 * 12 *48 ) » U(2*12.48)* PHI(12*48)* SIGMA ( 12 *48 ) * 

1 RCAP ( 12*48 ) * VISC ( 12 *48 ) ♦ RH0I12*48)* RH0U(12*48) 

CALL CRISIS (ZZ*AAA*AMPS* NNN) 

M = 1 
L = 2 

READ INPUT TAPE 5 * 100 *NFI LEI . NTAPE1 * NFILE* NTAPE* KMAX. KINC 

100 FORMAT (2IA/4I4) 

READ INPUT TAPE 5*101. EX* ZC. THETA 

101 FORMAT ( E10.3) 

THETAR = THETA*2. 0*3. 1416/360.0 
CALL LOCATE ( NFI LEI ♦ NTAPE1 ) 


READ 

TAPE NTAPE1* 

K* DI AM ( 2 ) * Z. AMPS. E 

(2) * W *QX 

% HWALL ( 2 ) 

1 

*PX 

» NMESH* P (2) * EZOjLOC»EPS 

*DWi DP* 

DZ 

READ 

TAPE NTAPE1.NI* i 

1 ( H <M * J . I ) . U (M * J* I ) ♦ RHOU 

( J.I). 

I = 1 • N I ) • 

1 

J=1 *NMESH) 




P ( 1 ) 

= P ( 2 ) - DP 




DO 11 

K = 3 *KMAX 




11 

HWALL (K) = 

HWALL ( 2 ) 




K = 2 

DIA * D I AM ( 2 ) 
CALL RHOAVE 


CALL STATEP 
QR = QX*PX/ 100.0 
Q = QX - QR 

CALL LOCATE (NFILE»NTAPE) 

CALL OUTPT 
DO 6 K = 3 *KMAX 

DZ = EX *DZ 
Z * Z + DZ 

IF(Z-ZC) 50*50.51 

50 DIAM(K) = D I AM ( 2 ) 

GO TO 52 

51 DIAM(K) s D I AM ( 2 ) + 2.0*(Z-ZC)*TANF(THETAR) 

52 CONTINUE 

DIA = DIAM(K) 

M = 1 
L = L + 1 
IFC3-L) 1.1*2 

1 M = 2 
L = 1 

2 CONTINUE 



CALL ENERGY 
IFIZ-ZC) 60*60*61 

60 CALL ITER 
GO TO 62 

61 CALL ITERS 

62 CONTINUE 

IF(NCHOKE) 4,4.3 

3 WRITE OUTPUT TAPE 6. 202* K* DW, U(M,2*2) 

GO TO 8 

4 CONTINUE 
CALL STATEP 
NN = NN - 1 

I FINN) 5,5,6 

5 NN a KINC 
CALL OUTPT 

6 CONTINUE 
8 CONTINUE 

REWIND NTAPE1 
REWIND NTAPE 
REWIND 8 

202 FORMAT ( 1H0, 18HFLOW CHOKED AT K = * 14, 10X* 

1 20HFLOW RATE ERROR IS • 2PF12.7* 10H PERCENT 

2 14 HCL VELOCITY = .0PF10.1, 8H FT/SEC ) 

CALL EXIT 

END 



VO 

Ov 


CHT0781 VAL WATSON SUBROUTINE ITERS FOR HT0780 
SUBROUTINE ITERS 

COMMON AMPS* W* WW. DP. DZ» EPS. FZO. Z. Q» QR * DW. DIA. EX. EXX. 
1 H. PHI. SIGMA. RCAP ♦ RHO. VISCt HWALL. RHOAV. RHOU. PHIW. 

1 P. E. DIAM. U. 

1 NMESH . K . L* M. LOC. NFILE. NTAPE . NCHOKE. NX. NNN 

DIMENSION 0 I AMI 100 1 * P(100l. E(100|. HWALL ( 100 ) 

DIMENSION HI2.12.48). Ut2.12.48). PHK12.48). S IGMA ( 12 .48 ) . 

1 RCAP (12.48) . VISCt 12 .48 ) . RH0(12»48)» RHOU(12.48) 

1ND = 0 

CALL CRISIS(ZZ.AAA) 

NCHOKE = 0 



DDP 

a DP 





P ( K ) 

= P 

7C 

1 

*-• 

+ 

DP 



LOC 

= 0 





NNN 

= 0 





DO 4 

N = 

1*30 





LOC 

= N 





CALL MOM 




NNN 

= 1 






CALL WDOT 





DW 

= tww - 

W)/W 




I F ( 

ABSFtDW 

) - EPS) 

10*20*20 

20 


I F ( 

I ND ) 22 

*22*1 


1 



DDP = 

DDP/2.0 


22 


I F ( 

DW) 3,10.2 


2 



DP = 

DP + DDP 





P(K) 

- P(K-l) 

+ DP 




GO TO 

4 


3 



DP = 

DP - DDP 





P(K) 

= P(K-l) 

+ DP 


IND = 1 


4 CONTINUE 
NCHOKE = 1 
10 CONTINUE 
RETURN 
END 



APPENDIX D 


FORTRAN PROGRAM TO EVALUATE GAS PROPERTIES USING PREPARED TABLES STORED 

ON A MAGNETIC TAPE 


These programs move prepared gas tables from a magnetic tape to the core 
storage and make interpolations (either logarithmic or linear) from the tables 
to obtain the gas density, thermal conductivity potential, electrical conduc- 
tivity, radiance, and viscosity from known values of enthalpy and pressure. 

The input for these programs is the magnetic tape prepared by the 
program in appendix E. 
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VO 

Oo 


CHT0767 VAL WATSON SUBROUTINE FOR AIR TABLES - INPUT ENTHALPY AND P 
SUBROUTINE NTAB (P, H, PHI* SIGMA* RCAP* VISC* NTAPE * NERR) 
DIMENSION PHIT ( 2 *300 ) * SIGMAT ( 2 *300 ) * RCAPT ( 2 * 300 ) * VI SCT (2*300) 
CALL CRISIS (ZZ *PHI T ) 

IF (IN) 1*1*2 
1 CALL LOCATE ( 2 * N T A P E ) 

READ TAPE NTAPE* (( SIGMAT ( I rJ ) * J = 1 * 281 ) * 1 = 1*2)* 

1 (( PHIT ( I *J) * J = 1 * 281 ) ♦ 1 = 1*2) * 

1 (( RCAPT ( I * J ) * J = 1 * 281 ) » 1 = 1,2) * 

1 (( V I SCT ( I ♦ J ) * J = 1 *281 ) * 1 = 1*2) 


REWIND NTAPE 
IN = 1 

2 NERR = 0 

IF (2.0E9-H ) 3*4,4 

3 NERR = 1 
GO TO 10 

4 IF (H - 2 *0 E7 ) 5*5*60 

5 I = XINTF (H/2»0E+5 ) + 1 

HH = (H/2.0E+5) - FLOATF(I-l) 

GO TO 7 

60 IF (H-2,0E8) 61,61,62 

61 I = XINTF (H/2 *0E+6 ) +91 

HH = (H/2.0E+6) - FLOATFU- 91) 

GO TO 7 

62 I = XINTFIH/2.0E+7) +181 

HH = (H/2.0E+7) - FLOATF( 1-181) 

7 PI = PH IT f 1 * I ) + ( PHI T( 1 ♦ 1 + 1 ) - PHI T ( 1 * I ) ) #HH 

P2 = PHIT (2*1) + ( PHITt 2*1 + 1 ) - PHI T ( 2 * I ) ) *HH 

51 = SIGMAT (1*1) + (SIGMAT (1*1 + 1) - SIGMAT ( 1 * I) I #HH 

52 = SIGMAT (2*1) + (SIGMATf 2 * 1 + 1 ) - SIGMAT (2*1)) #HH 

R1 = RCAPT (1*1) + ( RCAPT ( 1*1+1) - RCAPT (1*11) *HH 

R2 = RCAPT (2*1) + ( RCAPT (2*1+1) - RCAPT (2*1)) *HH 

VI = VI SCT (1,1) + ( V I SCT (1*1 + 1) - VI SCT (1*1)) *HH 

V2 = VI SCT (2*1) + ( V I SCT (2*1 + 1) - VI SCT ( 2 * I ) ) *HH 

PP = LOGF ( P ) / 2 • 3026 

PHI = PI + (P 2 - PI ) *PP 
SIGMA = EXPFtSl + ( S2 - S1)*PP) 

RCAP = EXPF(R1 + ( R2 - R1)#PP) 

VISC = VI + ( V2 - V1)*PP 
10 RETURN 
END 



CHT0768 VAL WATSON SUBROUTINE FOR AIR DENSITY TABLE 
SUBROUTINE NRHO (P* H* RHO* NTAPE* NERR) 

DIMENSION RHOT ( 2 #300 ) 

CALL CRISIS (ZZ*RHOT) 

IF (IN) 1*1*2 

1 CALL LOCATE (1 »NTAPE ) 

READ TAPE NTAPE. (( RHOT ( I * J ) * J-l *281 ) • 1*1*2) 
REWIND NTAPE 
IN = 1 

2 NERR = 0 

IF (2.0E9-H ) 3*4*4 

3 NERR = 1 
GO TO 10 

4 IF (H - 2.0E7) 5*5*60 

5 I * X I NTF (H/2*0E+5 ) + 1 

HH = (H/2.0E+5) - FLOATF(t-l) 

GO TO 7 

60 IF (H-2.0E8 ) 61*61*62 

61 I a XINTF(H/2,0E+6) +91 

HH = (H/2*0E+6) - FLOATF ( I- 91) 

GO TO 7 

62 I » XINTF(H/2*0E+7) +181 

HH = (H/2*0E+7) - FLOATF ( 1-181 ) 

7 RH1 a RHOT (1*1) + I RHOT (1*1+1) - RHOT(l*I))*HH 

RH2 a RHOT (2*1) + ( RHOT (2*1+1) - RHOT (2*1)) #HH 

PP a LOGF ( P ) /2* 3026 
RHO = P/(RH1 + ( RH2 - RH1)*PP) 

10 RETURN 
END 
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CHT0769 VAL WATSON SUBROUTINE FOR AIR TABLES - INPUT ENTHALPY AND P 
SUBROUTINE NTEMP ( P »H *SIGMA *NTAPE *NERR ) 

DIMENSION SI GMAT (2*300) 

CALL CRISIS ( ZZ »S IGMAT ) 

IF (IN) 1 > 1 • 2 

1 CALL LOCATE ( 3 iNTAPE ) 

READ TAPE NTAPE » ( ( SIGMAT ( I * J ) * J=l*281)* I»l*2) 

REWIND NTAPE 
IN b 1 

2 NERR * 0 

IF (2.0E9-H ) 3*4*4 

3 NERR = 1 
GO TO 10 

4 IF (H - 2.0E7) 5*5*60 

5 I = X I NTF (H /2.0E+5 ) + 1 

HH = (H/2.0E+5) - FLOATF(I-l) 

GO TO 7 

60 IF (H-2.0E8) 61*61*62 

61 I = X I NTF (H/2 *0E+6 ) +91 

HH = ( H/2 *0E+6 ) - FLOATF ( I- 91) 

GO TO 7 

62 I = X I NTF ( H/2 *0E+7 ) +181 

HH = (H/2.0E+7) - FLOATF (1-181) 

7 SI = SIGMAT (1*1) + ( S IGMAT (1*1 + 1) - S IGMAT ( 1 * I ) ) *HH 
S2 = SIGMAT (2*1) + ( S I GMAT ( 2 * I +1 ) - SIGMAT (2*1)) *HH 
PP=LOGF (P) / 2.3026 

SIGMA = (SI + (S2 - Si ) *PP ) 

10 RETURN 
END 



APPENDIX E 


FORTRAN PROGRAM FOR PREPARING GAS TABLES FOR USE IN THE PROGRAM 

IN APPENDIX D 


This program prepares the magnetic tape required for the program in 
Appendix D from any gas tables wherein the gas properties are given in terms 
of any two state properties* The program fits a third-order polynomial curve 
between the middle two points of each set of four adjacent points of the 
input tables. The gas properties are taken from this curve in equal incre- 
ments of enthalpy to form the table that is put on magnetic tape for use with 
the programs of appendix D. To check the fitted curves, the program plots 
the input data (shown as circles) and the fitted curves on the same graph, 
as shown in figure 2. (When the curve did not go through the symbols, the 
automatic plotter was out of adjustment; within the program the curve was 
forced through each point.) 
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CHT0766 PEGOT GENERATES TABLES AND PREPARES TAPE AND PLOTS 2401 

DIMENSION HI (100) *P ( 10) * PHI 1(100) • S I GMA1 ( 100 ) * RCAP1 ( 100 ) * V ISC1 ( 100H T0766 
1) »RH01(100) ♦HTEMPl(lOO) HT0766 

DIMENSION HA (300) *PHI A ( 300 ) *S IGMAA (300) *RCAPA< 300) * V ISCA ( 300 ) »RHOAHT0766 
1(300) *TEMPA (300) >YY(300 ) * RCAPLA ( 300 ) » RRHOA ( 300 ) HT0766 

DIMENSION H2( 100) * PH I 2 ( 1 00 ) »$IGMA2( 100) ♦ RCAP2 ( 100 ) * V ISC2 ( 100 ) *RH02HT0766 
1(100) >hTEMP 2< 100) *PTEMP1 ( 100 ) * PTEMP2 ( 100 ) ♦ STEMP1 ( 100 ) * STEMP2 ( 100 ) *HT0766 
2RTEMP1 (100) *RTEMP2( 100 ) *VTEMP1( 100) » VTEMP2 ( 100 ) tOTEMPl (100) *OTEMP2HT0766 
3(100) HT0766 

DIMENSION RHOT ( 2 *300 ) » TEMPT (2*300) » HT { 2 » 300 ) 

DIMENSION PH IT ( 2 * 300 ) » S I GMAT ( 2 » 300 ) » RCAPT ( 2 * 300 ) * VI SCT (2*300) 

CALL CRISIS ( 2Z » HI ) 

CALL LOCATE (1*8) 

P(l)®1.0 
P ( 2 ) e 10 .0 
P ( 3 ) = 0.5 
P ( 4 ) * 5.0 
P ( 5 ) = 50.0 
P ( 6 ) = 1.0 . 

C 

WRITE OUTPUT TAPE 6*299 
299 FORMAT ( 1H1 * 5X * 10H PAGE 1 ) 


1=1 HT0766 

CALL GENTAB ( HI .HTEMP1 *NXY1H .HA* TEMPA * DXX *N *0 * 2 *PHI A * TEMPA * P * I ) HT0766 

CALL GENTAB (STEMP1.SIGMA1 *NXY1S *TEMPA *S IGMAA *DXX»N*0*3*PH1A*TEMPA*HT0766 
1 P * I ) 

CALL GENTAB (RTEMP1.RCAP1. NXY1R* TEMPA *RCAPA* DXX *N *0*3 *PHIA* TEMPA *P*HT0766 
II) HT0766 

CALL GENTAB (VTEMP1*VISC1*NXY1V*TEMPA*VISCA*DXX*N *0*2 *PHIA*TEMPA*P*HT0766 
II) HT0766 

CALL GENTAB ( PTEMP 1 *PH 1 1 *NXY1P *TEMPA * YY *DXX *N * 1 *2 *PH I A* TEMPA *P * I ) HT0766 

CALL GENTAB (0TEMP1.RH0 1.NXY10* TEMPA *RHOA»DXX»N iO»4iPHI A »TEMPA»P* I (HT0766 
END FILE 7 

GO TO 300 HT0766 

301 WRITE OUTPUT TAPE 6»298 
298 FORMAT (1H1.5X. 10H PAGE 2 ) 

CALL GENTAB (H2 »HTEMP2 *NXY2H »HA» TEMPA »DXX*N*0i2*PHIA»TEMPA»P»I ) HT0766 

CALL GENTAB (STEMP2*SIGMA2 *NXY2S» TEMPA *S IGMAA* DXX *N*0»3» PHI A.TEMPA>HT0766 
IP • I > HT0766 

CALL GENTAB (RTEMP2 *RCAP2 *NXY2R ♦ TEMPA * RCAPA* DXX iN *0 *3 *PH I A • TEMPA *P»H T0766 
II) HT0766 

CALL GENTAB ( VTEMP2 *VI SC2 .NX Y2V * TEMPA* VISCA * DXX *N *0 * 2 *PH I A *TEMPA *P*H T0766 
II) HT0766 

CALL GENTAB(PTEMP2*PHI2*NXY2P*TEMPA.YY*DXX*N*1*2*PHIA*TEMPA*P.I ) HT0766 

CALL GENTAB (0TEMP2*RH02* NXY20* TEMPA *RH0A*DXX*N*0 *4*PHI A *TEMPA*P* I >HT0766 
END FILE 7 
300 DO 320 J* 1 * N 


HT0766 



noon no 


HT( I » J) =HA ( J ) 

TEMPT ( I » J ) =TEMPA(J) 
RHOT(IiJ) =RHOA ( J ) 
SIGMATt I * J ) =SIGMAA(J) 
RCAPT ( I * J) =RCAPA(J) 
VISCT ( I * J) sVISCA(J) 
320 PHIT ( I * J ) =PHIA( J) 


PREPARING VALUES FOR TAPE 

1=1+1 HT0766 

I F ( 1-2 ) 301 »301 *302 HT0766 

302 WRITE TAPE 8 ♦ (< RHOT ( I * J ) * J*1 »N ) • I =1 *2 ) HT0766 

END FILE 8 

WRITE TAPE 8 * ( ( SIGMAT ( I » J ) » J = l> N ). 1-1*2) * 

1 (( PHI T ( I » J ) * J = l* N )* 1 = 1*2) * 

1 {( RCAPT ( I » J ) * J=l* N )• 1=1*2) * 


1 ( ( VISCT ( I • J ) * J=1 • N )* 1=1*2) 

END FILE 8 

WRITE TAPE 8* ( ( TEMPT ( I . J ) ♦ J=1 * N ) *1 = 1*2) 

END FILE 8 
REWIND 8 

100 FORMAT ( 5E12 *5/ ( 5E12 • 5 ) ) 

PRINTOUT AND PLOTS 

CHECK FOR ERROR IF OVERRUN TABLES 
N£RR*NERR 
IF (NERR) 60*60*61 

61 WRITE OUTPUT TAPE 6.102 
GO TO 62 

60 WRITE OUTPUT TAPE 6* 103 

62 NERR = 0 
HOM = 2.5E6 
POM = 1.0 

CALL NRHO(POM.HOM*RHOA.8*NERR) HT0766 

IF (NERR) 70.70*71 

70 WRITE OUTPUT TAPE 6*102 
GO TO 72 

71 WRITE OUTPUT TAPE 6*103 

72 NERR = 0 

CALL NTAB ( POM*HOM*PHI A.SIGMAA *RCAPA *VISCA*8 *NERR ) HT0766 

IF (NERR) 80*80*81 

80 WRITE OUTPUT TAPE 6*102 
GO. TO 90 

81 WRITE OUTPUT TAPE 6*103 
90 CONTINUE 



H 

O 
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102 FORMAT (1H0* 30H ERROR SIGNAL INOP ) 

103 FORMAT ( 1H0 • 30H ERROR SIGNAL OKAY ) 

405 READ INPUT TAPE 5*400*JAM HT0766 

404 READ INPUT TAPE 5*400*NOC 

READ INPUT TAPE 5*2*XRANGE 
2 FORMAT (2E10.3) 

400 FORMAT ( 14) 

KATH = 1 
403 IN =0 

NPLOT = 0 
DO 401 K»2*6 

CALL NPLOTS(XNAME*YNAME*XRANGE*PfK) *CONVX *CONVY *SCALEX *SCALEY ♦ 

1 NPLOT * I N *ORIGX *ORIGY * CONVL ) 

I N — 1 

401 NPLOT=l 
407 KATHsKATH+1 

NOC»NOC-l 

IF ( NOC ) 402 *402 » 403 

402 END FILE 7 
JAMaJAM-1 

REWIND 8 

I F ( JAM ) 405*405*404 


HT0766 
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O 

ON 


CHT0785 PEGOT SUBROUTINE TO GENERATE TABLE VALUES , 401 

SUBROUTINE GENTAB ( X * Y *NXY *XX * YY * DXX *N *NUMB» JUMB * PHI A *TEMPA *P* I ) 
DIMENSION X ( 300 > *Y(300) ♦ XX (300) *YY(300) *PHI A ( 300 ) *TEMPA ( 300 ) * P i 10 ) 

1 «HA ( 300 ) 

CALL CRISIS (ZZZ.X) 

READ INPUT TAPE 5 * 199 *NAX .NAY 

199 FORMAT t 2A6 > 

READ INPUT TAPE 5*200»NXY 
READ INPUT TAPE 5.201*CX,CY 
READ INPUT TAPE 5. 201. OX. OY 

301 READ INPUT TAPE 5,201* (X( JJ)*Y( JJ) , JJel *NXY ) 

200 FORMAT ( 14 ) 

201 FORMAT ( 2E10 .3 ) 

WRITE OUTPUT TAPE 6 * 198 *NAX ,NAY ,CX ,CY * ( X (JJ ) *Y ( J J) ♦ J J=1 *NXY ) 

198 FORMAT (1H1 , 25X . A6 . 15X , A6//24X ♦ 10HC0NVERSI0N * 10X • lOHCONVERSION/ 

1 17X»2E20.8//(17X*2E20.8) ) 

NAX »NAY * =NAME OF X AND Y ARRAYS 
CX*CY ,=CONVERSION FACTORS 
OX*OY ,=0RIGINS FOR X.Y PLOTS 
X=X ARRAY OF INPUT TABLES 
Y=Y ARRAY OF INPUT TABLES 

NXY=NUMBER OF ARRAY VALUES IN EACH X AND Y 
XX=ABSCI SSA TABLE VALUES < IN OUR CASE HT) 

YY=ORDINATE TABLE VALUES I TEMPT *PHI T *SIGMAT *ETC. ) 

DXX=DELTA HT VALUE 

N=NUMBER OF VALUES IN FINAL TABLE (191) 

JUMB FOR UNITS NEED CORRECTED FOR TAPE 
NUMB=1 FOR PHIT TO USE SIMP SUBROUTINE*=0 FOR ALL OTHERS 
DMONsO.O 

GO T0(304. 302*310*312) *JUMB 
310 DO 331 JJal.NXY 
Y ( JJ ) *Y ( JJ ) *CY 
331 X( J J ) =X { JJ)*CX 

CALL SPLOTS ( X* Y*NXY*7.0*10.0*OX*OY*SX*SY*17) 

AY=LOGF ( CY) 

DO 313 JJ=1*NXY 
X(JJ)=X(JJ) /CX 

313 Y ( JJ )*LOGF<Y( JJ)/CY) +AY 
GO TO 304 

312 DO 314 JJ=1*NXY 

314 Y { JJ >»PU)/(Y( JJ) *CY) 

GO TO 304 

302 DO 303 JJ=1*NXY 
Y ( JJ )=Y( JJ >*CY 

303 X( JJ )=X( JJ )*CX 

304 YY ( 1 ) =Y ( 1) 

DXX*2.0E5 
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N=281 

XX(1)=0.0 
DO 305 J = 2,N 
IF (KITE ) 350.349.350 

349 IF (XX(J-l) - 1.99E7 ) 315,316.316 

316 DXX=2.0E6 

IF (XX(J-l) - 1.99E8 ) 315,317,317 

317 DXX=2.0E7 

315 XX ( J)=XX( J-ll+DXX 

350 CALL TAI NT( X »Y »XX I J),YY( J ) »NXY » 3 *NERR*DMON ) 

NERR=NERR 

GO TO (305,306,306) ,NERR 

306 WRITE OUTPUT TAPE 6,307,XX( J),YY( J) 

307 FORMAT ( 1H1 , 28X .25HERROR IN TAINT SUBROUTINE //30X»2HXX»12X,2HYY/24 
1X.2E15.8) 

CALL EXIT 
305 CONTINUE 
KITE = 1 

I F ( NUMB) 308,351,308 

308 PH I A ( 1 ) =YY ( 1) 

DO 311 J = 2 ♦ N 

CALL S IMP ( PH I A ( Jl.TFMPA , YY ,J,II) 

GO TO (311, 309, 309*309). II 

309 WRITE OUTPUT TAPE 6,300,1 1 , J 

300 FORMAT ( 1H1 , 28X , 25HERROR IN SIMP SUBROUTINE //30X.2HI .12X.2HJ /24 
IX, 2110) 

CALL EXIT 
311 CONTINUE 

351 GO TO (321,321,332,321) ,JUMB 

332 DO 333 K=1»N 

333 Y(K) = EXPF ( YY ( K ) ) 

CALL PLOTWS(OX,OY,SX,SY»XX. Y,N.7»-17»NERR ) 

GO TO 322 

321 CALL SPLOTS ( X , Y ,NXY ♦ 7. 0 , 10.0 ,OX ,OY , SX , SY , 17 ) 

CALL PLOTWS ( OX ,OY ,SX ,SY ,XX , YY ,N ,7 ,-17 ,NERR ) 

322 WRITE OUTPUT TAPE 6 ,253 ,NAX ,NAY ,SX ,SY ,OX ,OY 

253 FORMAT ( 1H1, 5X.A6, 5X ,A6 , 5X , 10HSCALE X = »E10.3»5X,10HSCALE Y = 

1 . E10 . 3 » 5X » 10H ORIG X = .E10.3.5X.10H ORIG Y = ,E10»3 ) 

RETURN 


A 
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401 


CHT2136 EVA PF60T SET UP SCALES AND PLOT 

SUBROUTINE SPLOTS(X*Y.NUM*RANGEX*RANGEY*ORIGX*ORIGY. SCALEX .SCALEY* 
1 NBYM ) 

DIMENSION X ( 500 ) * Y(500) 

CALL CRISIS (ZZZ. All) 

XM*X ( 1 ) 

DO 11 J= 1 *NUM 
IF (XM-X( J) ) 10*10*11 

10 XM= X ( J ) 

11 CONTINUE 
YM=Y ( 1 ) 

DO 13 J*1 *NUM 

IF I YM-Y( J) ) 12*12*13 

12 YMsY(J) 

13 CONTINUE 

OY = YM/RANGEY - 1.0E-30 
EXPQY = ( LOGF (QY)/2. 30258 ) 

IF (EXPOY ) 1*1*2 

1 EXPQY = EXPQY -1.0 

2 EXPQY = INTF ( EXPQY ) 

AQY = EXPF (LOGF (QY) - EXPQY * 2.30258 ) 

IF (AQY -2.0 ) 32*34*34 
32 SCALEY =2.0*10«0**EXPQY 

GO TO 41 

34 IF (AQY - 4.0 ) 36*36*38 

36 SCALEY = 4.0*10.0**EXPQY 
GO TO 41 

3e SCALEY = 10*0*10. 0**EXPQY 

41 QX = XM/RANGEX -1.0E-30 

EXPQX = (LOGF (QX )/ 2.30258 ) 

IF (EXPQX ) 3*3*4 

3 EXPQX = EXPQX -1.0 

4 EXPQX = INTF ( EXPQX ) 

AQX * EXPF (LOGF (QX) - EXPQX * 2.30258 ) 

IF (AQX - 2.0 ) 42*44.44 

42 SCALEX = 2.0*10.0**EXPQX 

GO TO 51 

44 IF ( AQX - 4.0 ) 46*46*48 
46 SCALEX = 4,0*10.0*#EXPQX 

GO TO 51 

48 SCALEX = 10. 0*10*0** EXPQX 

51 CALL PLOTWS (OR I GX .OR I GY .SCALEX* SCALEY *X*Y*NUM»7*NBYM*NERR) 
RETURN 
END 
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APPENDIX P 


FORTRAN PROGRAM FOR PLOTTING IN PICTORIAL FORM THE RESULTS 
OF THE PROGRAM IN APPENDIX B 


This program reads from magnetic tape the results of the program in 
appendix B and plots these results in pictorial form as shown in figures 4l 
through 44. 
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CHT073T VAL WATSON PLOTS OUTPUT FOR HT0720 PLOTS TEMP 401 

DIMENSION H ( 100 ) » R(100)t RHOU(IOO)* U< 100 ) • RUH(IOO)* RUU(IOO)* 

1 HH(IOOO), HR ( 1000 ) t HAVE(1000)» RUHA(IOOO)* HRAVE(IOOO). 

2 QR ( 1000) * QT(IOOO)# PRAD(IOOO)* E(lOOO)* P(1000)* 2 1 1000 ) 

• »TEMP ( 100 ) 

CALL CRI S IS { ZZ *H ) 

2 READ INPUT TAPE 5* 104* KSt KL* Kit NTAPE 

READ INPUT TAPE 5* 100* CRA* CRB* CHt CRU# CRUH 

READ INPUT TAPE 5* 100* SR* SH* SRUt SRUH* SU* SRUU* EX 

READ INPUT TAPE 5tlOO*STEMP 

READ INPUT TAPE 5* 104* NOTHC 

IF(NOTHC) 400*400*402 

400 READ INPUT TAPE 5* 100* A* CPOK 

402 CONTINUE 

100 FORMAT (E10* 3 ) 

104 FORMAT (14) 

IF (NT -NTAPE ) 40*42*40 

40 CALL LOCATE { 1 *NTAP£ ) 

KEND = 1 
NT = NTAPE 
42 KSKIP = KS-KEND 

IF(KSKIP) 700*702*700 
702 KSKIP = -0 
700 CONTINUE 

CALL SKIP(KSKIPtNTAPE) 

KEND = KL 

NUMM = KL - KS + 1 
KR = KI 

DO 50 K a 1 *NUMM 

READ TAPE NTAPE* M. DIAM* Z ( K) • AMPS* E( K ) * W* QT ( K) * HWALL * 

1 PRAD(K)* NMESH * P(K>* FZO* LOC* EPS* DW* DP* DZ* 

2 HAVE ( K ) * RUHA(K) * HRAVE(K) 

READ TAPE NTAPE. ( R ( J ) * H( J) * U( J)* RHOU(J)* J*1*NMESH) 

U(NMESHP) = 0.0 
RHOUt NMESHP) » 0.0 
R ( 1 ) * 0.0 
NMESHP * NMESH + 1 
R ( NMESHP ) * DI AM/2.0 
H( NMESHP ) = HWALL 
HH ( K ) = H ( 2 ) 

HR ( K ) = HAVE ( K ) /HH ( K ) 

QR ( K ) = QT ( K )*PRAD ( K ) /100. 0 
KR = KR + 1 
IF(KI-KRI 10*10*12 
KR = 0 

DO 8 J=1*NMESHP 

RUU ( J J = RHOU ( J )*U ( J } 
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RUH(J) a RHOU(J)*H(J) 

8 CONTINUE 

DOR = EX*Z(K) 

ORA = CRA + DOR 
ORB = CRB + DOR 
OH = CH + 0.5*DOR 
ORU = CRU + 0.5*DOR 
ORUH = CRUH + 0«5*DOR 
NUM = 2*NMESH+1 
NP = NUM 

DO 20 JP = 1 »NMESHP 
JJP = (NMESHP-JP+1) 

R ( NP ) = R ( JJP ) 

H ( NP ) = H(JJP) 

U ( NP ) = U(JJP) 

RHOU(NP) = RHOU(JJP) 

RUU(NP) = RUU(JJP) 

RUHINP) = RUH(JJP) 

CALL MTEMP ( P ( K ) * H ( NP ) * TEMP ( NP ) * 8 *NERR ) 

NP » NP-1 

20 CONTINUE 

DO 22 JP=2*NMESHP 
JPN s (JP+NMESH) 

RINP) = O.O-R(JPN) 

H C NP ) = H(JPN) 

U ( NP ) » U(JPN) 

RHOU(NP) = RHOU(JPN) 

RUU(NP) = RUU(JPN) 

RUHINP) = RUH(JPN) 

TEMP(NP)aTEMP( JPN) 

NP a NP-1 

22 CONTINUE 

CALL PLOTWS ( ORA » OH * SR * SH . * R* H * NUM >7 *-l * NERR ) 

CALL PLOTWS (ORB * OH * SR • SU * Ri U *NUM*7 »-l*NERR) 

CALL PLOTWS ( ORA * ORU ♦ SR • SRU t R» RHOU » NUM *7 »-l #NERR ) 

CALL PLOTWS ( ORB » ORU ♦ SR » SRUU» R* RUU *NUM *7 *-l * NERR ) 

CALL PLOTWSfORA* ORUH* SR ♦ SRUH* R* RUH *NUM*7 *-l»NERR) 
CALL PLOTWS ( ORB * ORUH » SR *STEMP *R *TEMP *NUM*7 *-l * N ERR ) 

12 CONTINUE 

IF(K-NUMM) 24*50*50 
24 CALL SKI P ( 1 » NTAPE ) 

50 CONTINUE 

END FILE 7 

CALL SPLOTS(Z,hh* NUMM* 10,0* 7.0* -14. 0*-14. *SSX *SS1 *-l ) 

CALL SPLOTS(Z,HAVE* NUMM* 10.0* 7.0* -14. 0*-5 .0 *SSX *SS2 *-l ) 

CALL SPLOTS(Z,RUHA* NUMM* 10,0* 7.0* -14. 0*+4. *SSX *SS3 *-l ) 

CALL SPLOTSIZ.E* NUMM* 10,0* 7.0* 0. 0 *-14. *SSX *SS4*-1 ) 

H 

P 




112 




CALL SPLOTS(Z,QT* NUMM* 10.0* 7.0* 0.0*- 5. *SSX *SS5»-1 ) 

CALL PLOTWS ( 0. *-5.*SSX*SS5*Z*QR» NUMM *7 *-l*NERR) 

CALL SPLOTS(Z,HR* NUMM* 10.0* 7.0* 0.0* 4. .SSX *SS6 *-l ) 

CALL SPLOTS ( Z ♦ P * NUMM* 10.0* 7.0* 0.0* 4. *SSX*SS7*-1) 

IF(NOTHC) 500*500*502 
500 ZO = W*CPOK/3.1416 

ZBAR = Z ( NUMM ) / ( ZO ) 

HINF = 0.307 *CPOK*AMPS/(SQRTF( A)*DIAM/2«0 ) 

HINFA = HINF*0»433 

QINF = 0.3831 *AMPS/ ( SQRTF ( A ) #D I AM#D I AM/4.0 ) 

EINF = 2.40/ ( SQRTF (A)*DIAM/2.0) 

HOBAR = HH ( 1 ) /HI NF 
RUHAB = RUHA(l) /HINFA 
HAVEB = HAVE ( 1 ) /HINFA 
SCZ = SSX/tZO) 

SCI = SS1/HINF 
SC2 = SS2/HINFA 
SC3 = SS3/HINFA 
SC4 = SS4/EINF 
SC5 = SS5/QINF 

CALL TCURVE(ZBAR*RUHAB* -14.0* -14.0* SCZ* SCI* NERR ) 

CALL TCURVE(ZBAR*RUHAB* -14.0* -5.0* SCZ* SC2* NERR) 

CALL T CURVE ( ZBAR * RUHAB * -14.0* 4,0* SCZ* SC3 * NERR) 

CALL ECURVE (ZBAR * RUHAB * 0.0* -14.0* SCZ* SC4* NERR) 

CALL TCURVEtZBAR. RUHAB* 0.0* -5,0* SCZ* SC5 * NERR) 

502 CONTINUE 

END FILE 7 

CALL SIMP (VOLT *Z*E*NUMM*NERR) 

NERR=NERR 

GO TO (900*901*901*901) *NERR 

901 WRITE OUTPUT TAPE 6*902 *NERR 

902 FORMAT ( 1H1 * 25H ERROR IN SIMP NERR a *12 ) 

CALL EXIT 

900 WRITE OUTPUT TAPE 6* 204 

204 FORMAT ( 1H1 * 30HPLOTS IN MKS UNITS / 

11X* 50H SHEET 1 - RADIAL PROFILES / 

4 5X • 50H ENTHALPY - LOWER LEFT / 

4 5X. 50H MASS FLUX - CENTER LEFT / 

4 5X • 50H ENERGY FLUX - UPPER LEFT / 

4 5X • 50H VELOCITY - LOWER RIGHT / 

4 5X • 5 OH MOMENTUM FLUX - CENTER RIGHT / 

4 5X* 50H TEMPERATURE - UPPER RIGHT ) 

SL = 1.0/EX 
SR = ABSF(SR) 

WRITE OUTPUT TAPE 6* 210, SR, SL, SH* SRU, SRUH* SU* SRUU,STEMP 
210 FORMAT ( lHO* 40H SCALES FOR RADIAL PROFILES ARE / 

1 20X * 30HRADIAL DISTANCE 


, E10.3 / 



p 


1 20X * 30HAXIAL DISTANCE 
1 20X * 30HENTHALPY 
1 20X » 30HMASS FLUX 
1 20X» 30HENERGY FLUX 
1 20X ♦ 30HVELOCITY 
1 20X » 30HM0MENTUM FLUX 
1 20X • 30HTEMPERATURE 
WRITE OUTPUT TAPE 6* 200 
200 FORMAT ( 1 HO • 


1 50H SHEET 2 - 

2 5X * 50H 

3 5X« 50H 

4 5X * 50H 
4 5X» 50H 
4 5X * 50H 
4 5X * 50H 
4 5X» 50H 

WRITE OUTPUT TAPE 
• *VOLT 


AXIAL PROFILES 
CENTERLINE ENTHALPY 
SPACE AVERAGE ENTH 
MASS AVERAGE ENTH 
VOLTAGE GRADIENT 
WALL HEAT FLUX 
PRESSURE GRADIENT 
RATIO OF AVE TO CL 


LOWER LEFT 

- CENTER LEFT 

- UPPER LEFT 

- LOWER RIGHT 

- CENTER RIGHT 

- UPPER RIGHT 
ENTHALPY - UR 


E10.3 

E10.3 

E10.3 

E10.3 

E10.3 

E10.3 

E10.3 


6* 202* SSX, SSI* SS2* SS3 * SS4* SS5 * S S7> 


202 FORMAT ( 1H0* 40H SCALES FOR AXIAL PROFILFS ARE 
1 20X * 30HAXIAL DISTANCE 
1 20X * 30HCENTERLINE ENTHALPY 
1 20X. 30HSPACE AVERAGE ENTHALPY 
1 20X. 30HMASS AVERAGE ENTHALPY 
1 20X * 30HVOLTAGE GRADIENT 
1 20X * 30HWALL HEAT FLUXES 
1 20X* 30HPRESSURE GRADIENT 
1 20X * 30HRATIO OF AVE TO CL ENTHALPY 
1 20X ♦ 30HCONSTRICTOR VOLTAGE 
RADIUS=DIAM/2.0 

WOVA=W/< 3. 14159*DIAM*DI AM/4.0) 

HWOVA=HtNF*WOVA 

SHINF=SQRTF(HINF) 


* E10.3 
» E10.3 
l El 0* 3 
» E10.3 

* E10.3 

* E10»3 

* E10.3 
» E10.3 
» E10.3 


/ 

/ 

/ 

/ 

/ 

/ 

) 


/ 

/ 

/ 

/ 

/ 

/ 


) 

SS6 


/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

) 


PMKS»1.013E5*P(1) 

WRITE OUTPUT TAPE 6*3O0*ZO»RADIUS »HINF*WOVA»HWOVA»SHINF*PMKS* 

.pm 


300 FORMAT I1H1*4X»20H(LC)Z(0) 

. 5X.20HR 

. 5X * 20H ( LC ) H ( INF! 

. 5X *20H ( LC ) W/A 

. 5X,20H(LC)H( INF) (LCIW/A 

. 5X.*20H(LC)H( INF) ##1/2 

. . 5X *20H ( LC ) P ( 0 ) 

»5X*E12«5*15H ATM 
WRITE OUTPUT TAPE6.301 *ZO*RADIUS 

301 FORMAT! 1 HO *4X*20H(LC)Z(0) 

» 5X.20HR 


♦E12.5*15H M / 

*E12.5 * 1 5H M / 

♦ E12.5 *15H J/KG / 

»E12.5»15H KG/SM**2 / 

• E 12# 5 » 15H W/M*#2 / 

*E 12. 5 * 1 5H M/S / 

♦E12.5.15H NEWTONS/M** 2 > 

//) 

»HINF*HINFA»EINF*QINF»PMKS*P(1) 
*E12.5 » 15H M / 

.E12.5.15H M / 



TTC 


• 5X,20H(LC)H( INF) »E12«5*15H J/KG / 

. 5X*20HH( INF) >E12.5»15H J/KG / 

. 5X»20HE( INF) fE12.5*15H V/M / 

. 5X»20H(LC)Q( INF) *E12.5*15H W/M**2 / 

• 5X*20H(LC)P(0) »E12# 5 » 15H NEWT0NS/M**2 ♦ 

•5X*E12»5 *15H ATM //) 

WRITE OUTPUT TAPE 6»302 »CPOK#A 
302 FORMAT ( 1H0*4X*29H ( LC ) SIGNIFIES LOWER CASE 

• / 5X *29H( > SIGNIFIES SUBSCRIPT /// 

•5X»20HCP/K .E12.5*15H MS/KG / 

.5X.20HA »E12 »5 » 1 5H 1/V**2 ) 

GO TO 2 
END 


CHT0738 VAL WATSON - PLOT OF 1/F(Z) 

SUBROUTINE ECURVE ( ZBAR >HOBAR »OX »OY*SX *SY »NERR ) 
DIMENSION F ( 110 ) * Z(llO) * E f 110) 

NERR=0 

IF (l.O-HOBAR) 1 » 1 *2 

1 NERR = 1 
GO TO 10 

2 ZS=0.0- LOGF ( 1 • 0-H0BAR*H0BAR ) /1 1 . 5 
DZ = ZBAR/100.0 

ZR = ZS 
DO 3 1=1*101 

F ( I ) = SQRTFtl.O - EXPF(-11.5*ZR) ) 

Ed) = 1.0/F ( I ) 

IF(E(I)-5.0) 20 *20*22 
22 Ed) = 5.0 

20 CONTINUE 

Z(I) = ZR -ZS 

3 ZR = ZR + DZ 

CALL PLOTWS (OX*OY.SX.SY»Z*E*101*7*-17*NERR) 

10 RETURN 
END 


A 





H CHT0743 VAL WATSON - PLOT OF F(Z> 

^ SUBROUTINE T CURVE ( ZBAR *HOBAR *OX *OY *SX *SY *NERR ) 

DIMENSION F(llO) * ZI110) 

NERR=0 

IF (l.O-HOBAR) 1*1*2 

1 NERR = 1 
GO TO 10 

2 ZS=0.0- LOGF ( 1 • 0-HOBAR*HOBAR ) /1 1 « 5 
DZ = ZBAR/100.0 

ZR = ZS 
DO 3 1=1*101 

F ( I ) = SQRTFd.O - EXPF (-11* 5#ZR ) ) 

Z ( I ) = ZR -ZS 

3 ZR = ZR + DZ 

CALL PLOTWS(OX*OY*SX»SY»Z*F*101*7*-17*NERR) 

10 RETURN 
END 



REFERENCES 


1. Stine, Howard A.; and Watson, Velvin R. : The Theoretical Enthalpy 

Distribution of Air in Steady Flow Along the Axis of a Direct -Current 
Electric Arc. NASA TN D-1331, 1962. 

2. Marlotte, G. L. ; Harder, R. L. ; and Prichard, R. W. : The Radiating Arc 

Column. AGARDograph 84, pt . 2, 1964, pp. 633-672. 

3. John, R. R., et al. : Thirty -Kilowatt Plasmajet Rocket-Engine Develop- 

ment - Third-Year Development Program. Avco Rep. RAD-SR -64-80, 

March 1964- 

4. Weber, H. E. : Growth of an Arc Column in Flow and Pressure Fields* 

AGARDograph 84, pt. 2, 1964, pp. 845-881. 

5. John, R. R., et al. : Theoretical and Experimental Investigation of Arc 

Plasma -Generator Technology. Part I. Applied Research on Direct and 
Alternating Current Electric Arc Plasma Generators. Avco Rep. 
ASD-TDR-62-729, Sept. 1963. 

6. Eckert, E. R. G.; and Anderson, J. E. : Performance Characteristics of a 

Fully -Developed Constricted Transpiration-Cooled Arc. AGARDograph 

84, pt. 2, 1964, pp. 751-795- 

7. Watson, Velvin R. : Comparison of Detailed Numerical Solutions With 

Simplified Theories for the Characteristics of the Constricted-Arc 
Plasma Generator. Proc. 1965 Heat Transfer and Fluid Mechanics 
Institute, Stanford Univ. Press, 1965, PP- 24-41. 

8. Masser, P. S. : Arc Jet Design. ARS Paper 2352-62, March 1962. 

9. Ahtye, Warren F. ; and Peng, Tzy -Cheng: Approximations for the Thermo- 

dynamic and Transport Properties of High Temperature Nitrogen With 
Shock -Tube Applications. NASA TN D-1303, 1962. 

10. Yos, Jerrold M. : Transport Properties of Nitrogen, Hydrogen, Oxygen, 

and Air up to 30,000° K. Research and Advanced Development Division, 
Avco Rep. RAD TM 63-7, March 1963* 

11. Spitzer, L. : Physics of Fully Ionized Gases. Interscience Pub., Inc., 

N. Y., 1956. 

12. Nardone, M. C.j Breene, R. G.j Zeldin, S. S.j and Riethof , T. R. : 

Radiance of Species in High Temperature Air. Space Sciences Laboratory, 
General Electric Co., R63SD3, June 1963* 

13. Anderson, J. E. : Transpiration Cooling of a Constricted Electric-Arc 

Heater. Rep. AFARL 66-0157, Aerospace Res. Lab., Wright Patterson AFB, 
Ohio (Heat Transfer Lab., Univ. of Minnesota, AF contract no. 

AF 33(657) -7380) , March 1966. 


117 


TABLE I . - SOURCES OF THERMODYNAMIC AND TRANSPORT PROPERTIES 

FOR HYDROGEN AND NITROGEN 
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TABLE II . - SEVERAL INTERESTING PROPERTIES OF THE CONSTRICTED ARCS THAT ARE SHOWN PICTORIALLY IN FIGURES 





LEGENDS FOR THE COLUMNS OF TABLE II 


Figure number - figure in which the numerical solutions are shown pictorially 
at the end of this report. 

Gas - type of gas used in the calculations (the gas properties used are shown 
in fig. 3) . 

Constrictor diameter - diameter of the constant area portion of the 
constrictor. 

Current - the total arc current passed through the constrictor. 

Flow rate at exit - the total gas flow passed through the constrictor (when 
no gas enters through the constrictor walls , this is equal to the gas 
flow rate at the inlet). 

Transpiration cooling mass flux - the gas flow transpired through the 
constrictor walls. 

Inlet pressure - pressure at the constrictor inlet. 

Constrictor length - the length between the constrictor inlet and the 
constrictor exit. 


Percent radiation at exit - the percentage of the wall heat flux at the 
constrictor exit that is due to the radiation heat losses. 


Mass average enthalpy at exit - the energy flux at the constrictor exit 

divided by the flow rate (this is the average enthalpy obtained by sub- 
tracting the total heat losses from the total power input and dividing 
by the flow rate - normally called the first law enthalpy measurement). 


Space average enthalpy at exit 


1 

A 


j> ^ 


at the constrictor exit. 


Voltage gradient at exit - the voltage gradient within the arc measured at 
the constrictor exit. 


Wall heat transfer rate at exit - the heat transfer rate impinging on the 
constrictor wall at the exit. 

Maximum wall heat transfer rate - the maximum heat transfer rate impingement 
on the constrictor wall. 

Length at maximum wall heat transfer rate - the axial position at which the 
maximum wall heat transfer rate occurs. 


120 



Computing time - the computing time required on the IBM 7094 for this 
numerical calculation. 

Flow choked at exit - a check in this column indicates that the flow is 
aero dynamically choked at the constrictor exit. 
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TABLE III . - THEORETICAL ESTIMATES OF THE- PROPERTIES OF THE PLASMA PRODUCED BY 
A CONSTRICTED-ARC PLASMA GENERATOR FOR ONE SET OF OPERATING CONDITIONS 


Plasma generator size 


Diameter of the constant area constrictor 0.00635 m 

Length of the constant area constrictor 0.l6m 

Diameter of the nozzle exit 0. 0318 m 

Operating conditions 

Current 500 A 

Arc voltage 620 V 

Mass flow rate (nitrogen) 0.0022 7 kg/s 

Chamber pressure 4.05x10 s N/m 2 

(4.0 atm) 

Resulting maximum constrictor wall heat transfer rate . . . 4.8X10 7 w/m 2 

Efficiency of constricted arc 53 percent 


Plasma properties at the constrictor exit 

Center -line enthalpy 

Center -line energy flux density . . . . 

Center -line velocity 

Center -line temperature 

I’ressure 

Center -line electron number density . . 
Center-line degree of ionization . . . 
Center -line electrical conductivity . . 

Plasma properties at the nozzle exit 


Size of uniform stream (less than 10-percent variation) 0.01 m 

Center -line total enthalpy 1.5xl0 8 j/kg 

Center -line energy flux density 3-5X10 8 w/m 2 

Center -line velocity 9X10 3 m/s 

I Center -line density . . 1.2X10" 4 kg/m 3 

Center-line electrical conductivity (order of 

magnitude only) 4xl0 3 l/fi-m 

Center -line Reynolds number per cm (order of 

magnitude only) 100 

Center -line magnetic Reynolds number per cm (order of 

magnitude only) 0.5 

Plasma properties in the stagnation region ahead of a blunt test body 

Enthalpy 1. 5x10 s j/kg 

Temperature . . ..... 14,500° K 

Pressure 2xl0 4 N/m 2 

(0.2 atm) 

Electron number density ....... 4xl0 ls electrons/cc 

Degree of ionization 75 percent 

Electrical conductivity 5X10 3 l/fi-m 


.... 1.4X10 8 j/kg 

.... 8 . 7 XIO 9 W/m 2 

.... 4.2X10 3 m/s 

.... 16, 000° K 

.... 2.2X10 5 N/m 2 

3.5X10 17 electrons/cm 3 
.... 55 percent 

.... 7XIO 3 l/fl-m 
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Cathode 


Constrictor 



Divergent nozzle 



A 



Temperature, *K Temperature, °K Temperature, *K 


(a) Hydrogen, p = 1 atm, (b) Hydrogen, p = 1 atm, (c) Hydrogen, p = 1 atm, elec 

enthalpy. pres sure/ density. trical conductivity. 


Figure 2.- Values of gas properties used to prepare tapes for the numerical calculations. 







Temperature, °K Temperature, °K Temperature, °K 

(g) Hydrogen , p = 10 atm, (h) Hydrogen, p = 10 atm, (i) Hydrogen, p = 10 atm, elec 

enthalpy. pressure/density. trical conductivity. 


Figure 2.- Continued. 



Temperature, °K 


Temperature, ®K 


Temperature, °K 


(j) Hydrogen, p = 10 atm, (k) Hydrogen, p = 10 atm, (Z) Hydrogen, p = 10 atm 

thermal conductivity. radiation. viscosity. 


Figure 2.- Continued. 
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T*mp«rQture, °K 

(p) Nitrogen, p = 1 atm, 
thermal conductivity. 


Temp€ratur«, °K 


(q) Nitrogen, p = 1 atm, 
radiation. 


T#mp#rofur*, °K 


(r) Nitrogen, p = 
viscosity 


Figure 2. - Continued. 
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(s) Nitrogen, p = 10 atm, 
enthalpy. 


(t) Nitrogen, p = 10 atm, 
pressure/density. 


(u) Nitrogen, p = 10 atm, elec- 
trical conductivity. 


Figure 2.- Continued. 



Thermal conductivity. 
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Nitrogen 
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Nitrogen 
p MO otm 
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(v) Nitrogen, p = 10 atm, 
thermal conductivity. 


(w) Nitrogen, p = 10 atm, 
radiation. 


(x) Nitrogen, p = 1 atm, 
viscosity. 
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Figure 2.- Continued. 
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(a) Enthalpy. 


0 


(d) Velocity. 




(c) Energy flux density. 


Figure 6.- Parabolic enthalpy and uniform velocity inlet distributions. 





p u z /2p 0 



I = 580 A 
m = .00353 kg/s 
R = .00635 m 
z 0 = 3.50 m 


h« = 1.89 

m/A = 27.1 

h^m/A- 5.2S 

*/K, = 1.37 

p 0 * 1.17 

H„= 8.17 

E„ = 819 

q„= 1.19 



(g) Space-average enthalpy. 



(h) Mass-average enthalpy. (k) Pressure. 


Figure 7.- Concluded. 
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(f) Center -line enthalpy. 
.8 - 



z/z 0 


(h) Mass -average enthalpy. 


Figure 8. 


P/Po 
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(i) Voltage gradient. 




0 .01 .02 .03 .04 .05 .06 

z/z 0 

(k) Pressure. 


Concluded. 




H s /H 



(f) Center-line enthalpy. 



(g) Space-average enthalpy. 



0 .01 .02 .03 .04 .05 .06 

z/z 0 


(h) Mass -average enthalpy. 


(i) Voltage gradient. 



z/z 0 

(k) Pressure. 


Figure 9* - Concluded. 
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1.2 


1.0 



(g) Space -average enthalpy. 



0 01 .02 .03 .04 

z/z 0 

(h) Mass-average enthalpy. 


Figure 10. 



(i) Voltage gradient. 



(j) Local heat transfer rate 



z/z 0 


(k) Pressure. 

Concluded. 
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(f) Center -line enthalpy. 



(g) Space-average enthalpy. 



0 .02 .04 .06 .08 .10 


z/z 0 

(h) Mass -average enthalpy. 

Figure 11. 


(i) Voltage gradient. 

1.2 - 

" Average heat flux 



(j) Local heat transfer rate 



z/z 0 

(k) Pressure. 

Concluded. 







_ ( I I I I I I I I I I 

(f) Center -line enthalpy. 
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(g) Space-average enthalpy. 



(h) Mass -average enthalpy. 


Figure 12. 
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(i) Voltage gradient. 



(j) Local heat transfer rate 



0 .01 .02 .03 .04 


z/z 0 

(k) Pressure. 


Concluded. 
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(a) Enthalpy. 




(b) Mass flow. 



-1.0 -.5 0 .5 l.0 v 
r/R 



(d) Velocity- 





-5 0 .5 1.0° 
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(e) Momentum flux. 


i = .00485 kg/s 
I = .00635 m 


h m * 1.54 x 10® J/kg 

m/A* 38.8 kg/sm 2 
t^m/A* 5.98x10® W/m 2 
VhZ* 1.24 xIO 4 m/s 
p 0 » 1.36 x 10® N/m 2 (1.35 atm) 

H«* 6.67 xIO 7 J/kg 

E«* 819 V/m 

q * 9.73x10® W/m 2 


(c) Energy flux density. 


Figure 13. - Comparison of the numerical solutions with experimental measurements; 1. 27 -cm- diameter 

constrictor; I = 473 A, m = 4.85 g/s. 



(f) Center -line enthalpy. 
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(g) Space-average enthalpy. 



0 .01 .02 .03 .04 .05 


z/z 0 

(h) Mass -average enthalpy. 


(i) Voltage gradient. 
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(j) Local heat transfer rate 



0 .01 .02 .03 .04 

z/z 0 


(k) Pressure. 


Figure 13 . - Concluded. 



/njh/h^m/A) pu/(rh/A 





159 



(f) Center -line enthalpy. 

.8 - 



(g) Space-average enthalpy. 
. 5 - 



0 .01 .02 .03 .04 .05 

z/z 0 


(h) Mass -average enthalpy. 


Figure l4 



(i) Voltage gradient. 


Average heat flux 
to second constric tor 

“ Average heat flux 



(j) Local heat transfer rate. 
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Figure 15. - Comparison of the numerical soluti 

constrictor; I 


0 (d) Velocity. 



h«,= 1.56 xlO B J/kg 

m/A = 7.58 kg/sm 2 

h 0O m/A= 1.18x10 s W/m 2 

Vtv = I.25XI0 4 m/s 

p 0 = 5.IIXI0 4 N/m 2 (.504 atm) 

H„o - 6.77 x|0 7 J/kg 

E„= 1636 V/m 

g„= I.97X10 7 W/m 2 


ms with experimental measurements; 0. 635 -cm- diameter 
= 2k0 A, m = 0.24 g/s. 


I = 240 A 
m = .000240 kg/s 
R = .00318 m 
z 0 = .238 m 



H q /H 



(f) Center -line enthalpy. 



(g) Space-average enthalpy. 
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(h) Mass -average enthalpy. 



(i) Voltage gradient. 



(j) Local heat transfer rate. 
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z/z 0 

(k) Pressure. 
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Figure 15*- Concluded. 



youh/h^ (rh/A) 


.2 

0 


(a) Enthalpy. 
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(cL) Velocity. 
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Figure l 6 . - Comparison of the numerical solutions with experimental measurements; 0 . 63 5 -cm -diameter 

constrictor; I = 210 A, m = 0.14 g/s. 
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(g) Space-average enthalpy. 
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(h) Mass -average enthalpy. 

Figure l6. 
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(a) Enthalpy. 




(c) Energy flux density. 


(d) Velocity. 



-1.0 0 o 

r/R 


(e) Momentum flux. 


I - 210 A 
m = .000369 kg/s 
R = .00318 m 
z 0 = .365 m 


h 00 = 1.37 x|0 8 J/kg 

m/A= 11.7 kg/sm 2 

h^m/A = I.59XI0 9 W/m 2 

VbZ = I.I7XI0 4 m/s 

p 0 = 5.72 x 1 0 4 N/m 2 C565 atm) 

Ha> » 5.92 xIO 7 J/kg 

E„,= 1636 V/m 

q„= l.73x I0 7 W/m 2 


Figure 17 . - Comparison of the numerical solutions with experimental measurements; 0. 635 -cm-diameter 

constrictor; I = 210 A, m = 0-37 g/s. 



(f) Center -line enthalpy. 




z/z 0 

(h) Mass-average enthalpy. 




(j) Local heat transfer rate. 



0 .05 .10 .15 .20 .25 .30 .35 .40 

z/z 0 


(k) Pressure. 



Figure 17 •- Concluded. 
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I = 183 A 
m = .000265 kg/s 
R = .00318 m 
20 = .262 m 


h„= 1.19x10 s J/kg 

m/A = 8.37 kg/sm 2 

h^m/A = 9.97x10 s W/m 2 

= I.09XI0 4 m/s 

p 0 = 4.76XI0 4 N/m 2 (.470 atm) 

Ha, = 5.16 xlO 7 J/kg 

E„= 1636 V/m 

I.5IXI0 7 W/m 2 


I 


Figure 18.- Comparison of the numerical solutions with experimental measurements; 0. 63 5 -cm -diameter 

constrictor; I = 183 A, m = 0.27 g/s. 



(f) Center -line enthalpy. 




(h) Mass -average enthalpy. 


Figure l8. 


(i) Voltage gradient. 




Concluded. 
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(c) Energy flux density. 

Figure 19«- Comparison of the numerical solutions with experimental measurements; 0. 635-cm-diameter 

constrictor; I = 180 k, m = 0.50 g/s. 
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(f) Center -line enthalpy. 



(g) Space-average enthalpy. 
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z/z 0 


(i) Voltage gradient. 



(j) Local heat transfer rate. 
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(k) Pressure. 


(h) Mass -average enthalpy. 


Figure 19.- Concluded. 
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(a) Enthalpy. 





(d) Velocity- 



-1.0 0 ° (e) Momentum flux. 
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I * 150 A 
rh = .000139 kg/s 
R=. 00318 m 
z 0 = .l38 m 
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Hoc = 

4.23XI0 7 

J/kg 

E«> = 
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Figure 20.- Comparison of the numerical solutions with experimental measurements; 0 . 635-cm-diameter 

constrictor; I = 150 A, m = 0.14 g/s. 
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Figure 21.- The constricted arc with no radiation losses (radiance set equal to zero). 
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(f) Center -line enthalpy. 
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(g) Space -average enthalpy. 
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(h) Mass -average enthalpy. 


Figure 21 
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(i) Voltage gradient. 



( j) Local heat transfer rate 
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z/z 0 

(k) Pressure. 


- Concluded. 
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(f) Center -line enthalpy. 
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(g) Space-average enthalpy. 



(h) Mass -average enthalpy. 


Figure 22. 
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Figure 23.- The constricted arc with radiation losses much greater than conduction losses. 
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Figure 24.- The numerical solutions for a 1.27-' 
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(cL) Velocity. 
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constricted arc with viscous forces included 
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(d) Velocity. 
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m * .000227 kg/s 
R = 00318 m 
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h„ * 3.25 x I0 8 J/kg 

m/A * 7.16 kg/sm 2 
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-cm constricted arc with viscous forces included. 



(f) Center -line enthalpy. 



(g) Space -average enthalpy. 
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(h) Mass -average enthalpy. 
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(k) Pressure. 


Figure 26.- Concluded. 
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(f) Center -line enthalpy. 



J. 

I 

! 


(g) Space-average enthalpy. 
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(h) Mass -average enthalpy. 


Figure 28 


(i) Voltage gradient. 
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Concluded. 
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(a) Enthalpy. 
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Figure 29*- The numerical solutions wi 
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(d) Velocity. 
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radial convection of energy and momentum. 
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(h) Mass -average enthalpy. 
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(k) Pressure. 


Figure 29- - Concluded. 
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Figure 32.- The nitrogen arc with transpiration-cooled constrictor walls. 
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(a) Constrictor inlet. (b) 0.0029 m from the constrictor inlet. 

Figure 35 •- The asymmetrical constricted arc with an axial flow of gas; enthalpy and mass flux as 

functions of radius and azimuthal position; I = 5^0 A; m = 0.00353 kg/s; R = 0.00635 m; z 0 = 3*50 m; 
hoo = I. 89 XIO 8 J/kg; m/A = 27.8 kg/sm 2 ; p 0 = 1.17X10 5 N/m 2 (l.l6 atm). 
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